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INTRODUCTION 


The  increasing  need  for  information  concerning  the  effects  of  aerosols  upon 
electro-magnetic  radiation  passing  through  media  containing  these  aerosols 
requires  continual  improvement  of  existing  processes.  As  our  understanding  of 
these  processes  improves,  the  need  for  sophisticated  numerical  codes  capable 
of  accurately  representing  such  theory  continually  increases..  It 'is  to  this 
end  that  computer  program  AGAUSX  (and  its  predecessors)  has  been  constructed. 

PROGRAM  AGAUSX 

Experience  acquired  over  a  period  of  several  years  with  various  versions  of 
the  single  scattering  codes  PGAUSS  and  AGAUS  has  revealed  that  users  of  such 
codes  often  tend  to  be  over  conservative  in  their  choice  of  the  number  of 
radii  at  which  Mie  calculations  are  performed.  Many  single  scattering  codes 
require  that  the  user  specify  some  input  parameter  such  as  the  total  number  of 
radii  to  be  used.  Unless  a  user  has  had  a  great  deal  of  experience  with 
choosing  such  parameters  for  various  types  of  size  distributions,  there  is  a 
tendency  to  use  many  more  radii  than  may  be  needed  to  obtain  results  at 
acceptable  levels  of  accuracy.  Since  the  overall  running  time  of  Mie  codes 
may  be  greatly  reduced  by  reducing  the  number  of  particle  sizes  treated,  it  is 
desirable  to  have  a  code  in  which  one  specifies  an  acceptable  level  of 
accuracy  and  which  then  uses  only  as  many  particle  radius  values  as  are  needed 
to  satisfy  that  requirement.  Rather  than  specifying  the  number  of  radii,  the 
user  of  AGAUSX  must  specify  a  "convergence"  level  DELTA.  The  quantity  DELTA 
represents  the  minimum  fractional  accuracy  that  is  acceptable  for  certain 
results  of  a  run.  In  operation,  AGAUSX  then  runs  through  the  Mie  calculations 
and  integrations  over  size  using  only  a  few  particle  radii  and  obtains  first 
estimates  of  extinction  coefficients,  etc.  The  size  intervals  are  then 
reduced  to  one-half  their  former  value  and  new  estimates  are  calculated.  If 
the  new  and  old  estimates  agree  to  within  (DELTA  x  100)  percent,  then  AGAUSX 
ceases  to  treat  additional  particle  sizes  and  proceeds  to  its  final 
calculations.  If  the  new  estimate  does  not  agree  with  the  "old"  one  to  within 
(DELTA  x  100)  percent,  the  size  interval  is  again  cut  in  half  and  the 
calculations  and  comparisons  are  repeated.  The  cycling  occurs  repetitively 


n 


until  either  the  convergence  criterion  is  satisfied  or  until  some  precoded 
maximum  number  of  values  (513  for  AGAUSX)  for  particle  radius  have  been  used. 

It  is  believed  that  AGAUSX  will  represent  a  substantial  advancement  in  aerosol 
modeling  because  it  contains  a  large  variety  of  size  distributions,  the 
possibility  of  modeling  mixtures  of  aerosols  with  different  optical 
properties,  treats  changes  in  particle  size  with  variable  relative  humidity, 
is  compatible  with  ACT"*'  and  will  probably  provide  substantially  shorter 
computation  times  than  most  existing  aerosol  codes.  In  addition  AGAUSX  has 
the  capability  of  producing  an  analytic  phase  function  when  quick  results  of 
intermediate  accuracy  are  desired. 

In  summary,  program  AGAUSX  is  a  versatile,  single  scattering,  state  of  the  art 
computer  code  capable  of  performing  the  following: 

(1)  Provisions  for  producing  (a)  phase  functions,  (b)  Legendre 
coefficients  or  (c)  scattering  fractions  compatible  with  ACT. 

(2)  Provisions  for  producing  an  analytic  phase  function  represented  by 
Legendre  coefficients. 

(3)  Provision  for  automatic  cycling  over  a  range  of  wavelengths  and 
averaging  of  results  over  such  a  range  of  wavelengths. 

(4)  Provision  for  treating  changes  in  hygroscopic  particle  sizes  which 
may  occur  with  variations  in  relative  humidity. 

(5)  Provision  for  treating  multicomponent  aerosols  differing  in  specific 
gravity  and/or  optical  properties. 

(6)  Availability  of  options  whereby  users  can  either  define  particle 
number  densities  or  let  them  be  calculated  from  mass  density  and  mass 
concentration. 


"^Atmospheric  Sciences  Laboratory  -  Chemical  Systems  Laboratory  -  TRADOC 
Systems  Analysis  Agency,  Battlefield/ Smoke  Obscuration  Code. 
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(7)  Availability  of  options  whereby  users  can  select  various  size 
distribution  models,  including: 

(a)  Log  Normal 

(b)  Double  Exponential 

(c)  Dermend jian' s  Model  C 

(d)  Power  Laws 

(e)  Khirgian  -  Mazin 

(f)  Modified  Gamma  whose  control  parameters  may  include  liquid 
water  content 

(g)  Four  Bimodal  Log-Normal 

(h)  Marshall  -  Palmer  rain  model 

(i)  User  Supplied 

(8)  Inclusion  of  an  internal  subroutine  to  provide  automatic  look  up 
and/or  interpolation  of  optical  constants  for  liquid  water  at  wavelengths 
between  0.35pm  and  200pm. 

(9)  Inclusion  of  built-in  internal  checks  to  warn  users  if  certain 
computed  quantities  seem  to  be  failing  to  converge  and  extended  error  messages 
associated  with  failures  in  the  Mie  routine. 

Every  effort  has  been  made  to  insure  that  program  AGAUSX  is  machine 
independent.  To  this  end  AGAUSX  has  been  written  in  ASCII  FORTAN,  and  is 
available  in  the  form  of  BCD  punched  deck.  Should  the  user  have  any 
questions,  discover  possible  inaccuracies,  or  simply  be  desirous  of  a  punched 
deck,  please  contact  R.  C.  Shirkey,  AUTOVON  258-4200  or  (505)  678-4200. 
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MIE  THEORY 


Mie  theory1  predicts  the  scattering  by  and  the  absorption  in  an  isolated, 
discrete,  homogeneous,  isotropic  sphere  of  diameter  D  with  a  known  complex 
refractive  index  n  =  m-ik  relative  to  the  surrounding  medium  and  illuminated 
by  monochromatic  radiant  energy  with  wavelength  X  in  the  surrounding  medium. 
The  theory  is  given  in  detail  in  standard  texts  and  need  not  be  repeated 
here.  Instead,  only  those  elements  of  theory  needed  for  an  understanding  of 
the  numerical  algorithms  used  in  AGAUSX  are  included. 

Scatterers  attenuate  beams  of  radiant  energy  by  scattering  some  of  the  energy 
into  directions  other  than  the  incident  or  forward  direction  and  by  absorbing 
some  of  the  incident  energy  within  the  body  of  the  particle.  The  combined 
effect  of  pure  scattering  by  the  particle  and  true  absorption  within  the 
particle  is  termed  extinction.  The  amount  of  extinction,  scattering  and 
absorption  by  a  single  particle  is  given  in  terms  of  corresponding  equivalent 
blocking  areas  or  cross  sections,  C  t,  Csca,  and  Cabg,  respectively.  These 
cross  sections  depend  only  on  the  refractive  index  of  the  particle  n  =  m-ik 
and  the  size  parameter  a  =  2irr/X,  where  r  is  the  particle  radius,  and  X  is  the 
wavelength. 

The  transmission,  T,  of  a  cloud  of  particles  of  geometric  depth,  d,  and  number 
density,  N,  is  given  by 


T 


(1) 


1G.  Mie,  1908,  " Cbnsiderations  on  the  Optics  of  Thrbid  Media,  Especially 

Colloidal  Metal  Sols,"  Ann  Phys,  25 
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with  the  optical  depth,  t,  given  by 


T 


Kext  * d » 


(2) 


where 


K. 


ext 


=  NC. 


'ext ' 


(3) 


The  balance  between  loss  by  scattering  and  loss  by  absorption  is  frequently 
characterized  by  the  albedo  of  single  scattering  0)o,  given  by 


oo 


0 


C  C 

sea  sea 


C  +  C  ,  C  „ 
sea  abs  ext 


(4) 


A  scatterer  with  u>Q  =  1  has  no  absorption  and  is  termed  a  conservative 
scatterer.  The  albedo  gives  the  probability  that  a  photon  encountering  the 
scatterer  will  be  scattered  into  some  direction  including  the  incident 
direction. 

Although  the  extinction  by  a  cloud  of  particles  is  correctly  given  by 
equations  (1)  and  (2),  two  implicit  assumptions  may  lead  to  improper  use  of 
the  equations.  The  optical  depth  t  in  equation  (2)  does  not  include  losses 
caused  by  absorption  in  the  medium  surrounding  the  particles.  This  assumption 
obviously  breaks  down  at  wavelengths  for  which  atmospheric  gases  absorb 
appreciably.  The  second  assumption  is  that  scattered  photons  never  return  to 
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the  incident  direction,  i.e.  ,  that  there  is  no  multiple  scattering.  This 
effect  becomes  increasingly  important  as  optical  depths  exceed  x  =  0.1  for 

broad  band  propagation.  » J  For  laser  wavelengths  the  effect  is  increasingly 

.  9  ^ 

important  as  x  goes  beyond  a  value  of  15.  ’ 

A  final  caution  should  be  noted  in  regard  to  absorption  within  the  particle. 
Although  absorption  within  the  particle  is  correctly  determined  by  the 
wavelength-dependent  imaginary  part  k  of  the  refractive  index  n,  the  explicit 
mechanism  which  causes  the  absorption  is  usually  not  specified.  Usually  the 
absorption  is  joule  heating  and  it  is  sometimes  necessary  to  account  for  the 
isotropic  black  body  radiation  emitted  by  the  scatterer  when  its  temperature 
rises  above  that  of  its  surroundings.  There  may  also  be  circumstances  when 
quantum  transitions  occur  in  the  scatterer  followed  by  emission  at  or  near  the 
same  wavelengths.  It  is  incumbent  on  the  user  of  the  numerical  algorithms 
presented  here  to  properly  include  these  effects  since  they  are  not 
automatically  accounted  for  in  these  algorithms. 

All  scattering  properties  of  spheres  are  computed  from  m  and  k,  and  through 
the  use  of  the  induced  electric  and  magnetic  multipole  moments  of  the  sphere 
an  and  bn,  respectively.  The  moments  are  given  by* 


'F' (na)'F  (a)-n  ¥  (na)'F'(a) 

„  -  n  n'-  J  J  r\K  J  (5) 

n  V'CnajC  (a)-n  ¥  (na)^'(a)’ 
n  n  n  n 


^S.  Hoijer,  1974,  Atmospheric  Attenuation  of  a  Light  Beam  Due  to  Scattering 
and  Absorption,  Research  Institute  of  National  Defense,  FOA-2-C-2659-E1- 
E2-E3-E4 

O 

0.  Steinvall,  1974,  Computed  MIE  Scattering  Properties  for  Laser  Wavelengths 
in  Various  Atmospheric  Madia,  Research  Institute  of  National  Defense,  F0A-2-C- 
2662-E1-E3 

*Note  that  n  is  used  as  a  subscript,  an  integer  index,  and  a  complex  index  of 
refraction  when  it  is  not  a  subscript. 
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and 


nr  (not)  ¥n  (a)  (na)  r  (a) 

bn  "  (na)  ^  (a}  (na)  ^  (a) 


(6) 


The  prime  denotes  differentiation  with  respect  to  the  argument.  The 
^(z)  and  5n(z)  functions  are  Ricatti-Bessel  functions  of  the  first  and  third 
kind,  respectively,  and  are  related  to  the  spherical  Bessel  functions  jn(z) 
and  nn(z)  by 


¥  (z) 
n  J 


z3nCZ), 


(7) 


and 


f 

K  (z)  =  zj  (z)  -  izn  (z)  =  Y  (z)  +  ix  O) , 
n  n  n  n  n  '  (.o; 


where 


j„(z)  =  &l/2  Jntl/2w, 


(9) 
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and 


c  /-Tr  ^  1/2  .. 

nn(z:i  *  ^  N 


n+1/2 


(z) 


(10) 


The  function  On+l/2^z^  the  half  integral  order  Bessel  function;  the 

function  Nn+i/2^z^  t^ie  half  integral  order  Neuman  function. 

The  extinction  cross  section  is  computed  from 


C 


ext 


xi 

2tt 


00 


I 

n=l 


(2n+l)  Re  (a  +b  ) , 
n  n 


(ID 


and  the  scattering  cross  section  from 


sea 


X2 

2tt 


I  (2n+l)  [|an|2+|bn|2]. 


n=l 


(12) 


The  various  cross  sections  are  the  basic  quantities  used  in  scattering 
problems,  but  they  are  not  the  quantities  usually  computed  directly  from  Mie 
algorithms.  Instead,  it  is  more  convenient  to  compute  dimensionless 

efficiency  factors  QgXt  and  Qsca>  which  depend  on  n,  k,  and  a,  and  which  are 
multiplied  by  the  geometrical  sphere  cross  section  to  obtain  the  true  cross 

O 

section  =  irr  Q^.  Thus, 
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2 


no 


Q 


ext 


I 

n=l 


(2n+l)  Re  (a  +b  ) 
n  nJ  ’ 


(13) 


and 


sea 


a  n=l 


l  (2n+1)  [|a  |2+|b^|2] . 


n 1 


(14) 


Although  the  cross  sections  account  for  the  energy  removed  from  the  forward 
beam,  they  do  not  give  any  information  about  where  the  scattered  photons  go. 
This  information  is  contained  in  scattering  amplitudes  and  intensity  factors 
which  relate  the  flux  density  scattered  through  an  angle  0  relative  to  the 
incident  flux  density.  There  are  two  amplitudes,  S^(9)  and  82(6 )>  and  the 
intensity  factors  i^(0)  and  ^(O),  which  correspond  to  light  respectively 
polarized  perpendicular  and  parallel  to  the  plane  of  the  scattering  defined  by 
the  direction  of  incidence  and  the  direction  of  scattering. 

The  intensity  factors  are  related  to  the  scattering  amplitudes  by 


ijO)  =  Is^e)!2, 


(15) 


19 


and 


i2(0)  =  |S2(0) 


The  amplitudes  come  from  the  multipole  moments  through 


s1  (0) 


Y  2n+l 

n(n+l) 
n=l  v  } 


[a  tt  (0)  +b  t  (0)1, 
n  n  n  n  J  J  ’ 


and 


s2(0) 


oo 


l 

n=l 


2n+l 

n(n+l) 


[Vn(e)*anVe>)* 


and  angular  factors  irn ( 0 )  and  xn(0)  defined  in  terms  of  associated 
functions: 


(16) 


(17) 


(18) 


Legendre 


^(9)  =  p„(cos  6)/ sin  6; 


(19) 


(20) 


T 

n 


(9) 


dP1 (cos  6) 
n  ' 

d0 


The  amplitudes  have  relative  phase  6  =  argS^  -  argS2» 
Alternative  expressions  frequently  used  are 


dP  (cos  0) 
V0)  =  d(cos  0)  ’ 


(21) 


and 


dTT  (0) 
n 

TT  (0)  =  cos0* tt  (0)  -  sin20  •  -r? - ,  (22) 

n  n  d(cos  0) 


where 


P  (cos  0) 
n 


i  ,n 

1  Cl  f  2  n  i  ^ 

-  (cos  0-1) 

on  i  j  nQ 

2  n!  dcos  0 


(23) 


These  functions  satisfy  the  following  recurrence  relations: 
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(24) 


V0)  .  cos  e  Vi(6)  . 


V2(e>> 


and 


Tn(0)  =  COS0  [TTn(0)-TTn_2]-(2n-l)sin20Trn_1(0)+Tn  2(0)  . 


(25) 


The  scattering  cross  section  measures  the  ability  of  a  particle  to  scatter 
light,  and  it  is  to  be  expected  that  Cgca  is  obtained  from  an  integral  over 
the  scattering  intensity  factors.  Equation  (12)  follows  from 


A2  1 

Csca  “  4tF  /  [ij(9)  +  i2(0)]  dcos0. 


(26) 


Although  the  intensity  factors  themselves  may  be  used  in  scattering 
calculations,  they  are  primarily  suited  for  computing  flux  densities,  and  it 
is  frequently  more  convenient  to  measure  and  compute  scattered  light  in  terms 
of  radiances.  Radiances  do  not  have  the  1/r^  dependence,  and  it  is  therefore 
unnecessary  to  know  the  distance  from  the  scatterer  to  the  detector  if  the 
detector  field  of  view  is  small  and  is  filled  by  the  scattering  cloud.  The 
phase  function,  p(0),  gives  a  radiance,  I,  scattered  into  the  0  direction  in 
terms  of  the  radiance  IQ  incident  on  the  particle.  The  phase  function  is 
dimensionless  and  is  defined  here  as 
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(27) 


PCS) 


2ttC 


ext 


[ijO) 


+  i2(0)]. 


The  normalized  phase  function  p(0)dft/4  gives  the  probability  of  a  photon 
being  scattered  through  an  angle  6  into  an  element  of  solid  angle 
dft  =  d<l>d(cos  0).  The  integral  of  the  normalized  phase  function  is  the  single 
scattering  albedo  u>0,  which  gives  the  probability  that  the  photon  is 
scattered  : 


~  2tt  1 

«0  =  4?  /  J  P(e)d<Mcose 

o  -1 


4ttC 


ext 


1 

/ 


-1 


[i1(0)+i2(0)]dcos6. 


(28) 


or 


co  =  C  /C 
o  sea  ext 


(29) 


The  phase  function  contains  a  sum  over  the  polarization  states  implicit  in  the 
i^  and  i2  intensity  factors,  and  is  thus  unsuitable  for  describing  the 
polarization  of  the  scattered  light. 

The  phase  function  can  also  be  represented  by  a  Legendre  series: 


n-1 

P(6)  =  l  5f  P  (cos0),  (30) 

l=o  *  * 
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where  the  Legendre  expansion  coefficients  are  given  by 


(2SL+1)  r 

2  /  P(0) P£(cos0)d(cos0) , 


(31) 


and  Pa(cos0)  are  the  usual  Legendre  polynomials. 

ANALYTIC  PHASE  FUNCTION 


A  new  analytic  phase  function  has  been  constructed  that  approximates  the 
computed  phase  function.  This  analytic  phase  function  is  comprised  of  two 
analytic  functions;  a  new  function  up  to  angle  0^  and  a  modified  Henyey  - 
Greenstein  analytic  function  for  angles  greater  than  0  This  new  analytic 
phase  function  is  called  the  Goedecke-Henyey-Greenstein  analytic  phase 
function  and  will  subsequently  be  referred  to  as  the  GHG  phase  function. 

The  GHG  phase  function  has  the  following  analytic  form: 


PIP) 


l-g‘ 


(l-2gy+g2)3/2 


+  (1-  2") 


(l-2gU+g2) 1/2 


], 


1  * 
1-g 


(l-2gp+g2) 3/2 


P>  Pr 
(32a, b) 

Pj*  \)>,  -1, 


where  y  =  cos0 ,  u>0  is  the  albedo  for  single  scattering,  and  P(y)  is  the  phase 
function  at  y.  The  bar  above  quantities  denotes  approximations.  Here,  g,  a, 
0,  y^  are  parameters  that  must  be  fixed  and  are  evaluated  as  follows: 
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Using  equation  32a  as  a  generating  function,  we  may  write  the  equivalent 
equation 


P(U)  =  %  l  (aN+l)gNP  (y),  (33a) 

N=0 


or 


00 

PM  =  l  VNtw)- 

N=0  1N 


(33b) 


where  the  P^(p)'s  are  the  Legendre  polynomials,  and  again  the  bar  above  the 
symbols  denotes  approximations.  Making  use  of  equations  32a  and  33a  we  now 
solve  for  a  and  3  by  forcing  the  equations  to  match  at  p  =  1.  This  is  done 
because  forward  scattering  is  usually  dominant  for  large  Mie  size 
parameters.  Solve  simultaneously 


PCD  =  [j  d+g)  +  (i-  f)(i-g)l» 

d-g) 


(34a) 
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(34b) 


co 


2 


(a  +  l)g2. 


In  equation  (34b)  we  have  forced  the  approximate  value,  to  equal  the  exact 
value,  o>2.  u>2  has  been  used,  as  opposed  to  as  it  was  found  that  the  match 
produced  by  =  a)Q(a+l)g  will  not  work  for  all  cases.  Since  we  expect  that 
for  p  <  p^  the  actual  phase  function,  and  therefore  the  HG  part  of  the  GHG 
phase  function,  is  small  compared  to  its  values  for  p  >  p^,  we  take  the  values 
of  a  and  g,  found  by  simultaneously  solving  equations  34a  and  b,  as  our 
working  values. 

Ideally  p^  could  be  found  by  setting  equation  34a  equal  to  zero;  however, 
using  a  value  of  p^  determined  in  such  a  manner  would  produce  an  unnatural  dip 
in  the  phase  function  where  we  match  32a  to  32b.  Therefore  we  take  p^  to  be 
that  value  of  p  for  which  the  phase  function  represented  by  equation  (32a) 
first  goes  negative  (the  phase  function  represented  by  equation  32a  may  not  go 
negative,  but  dependent  upon  the  various  parameters  it  may) .  3  is  found  by 
matching  32a  to  32b  at  p  =  p^.  This  yields 


3  -  |  +  (1-  (l-2gyi+g2)/(i-g2)  >  o. 


(35) 


Since  p^  is  now  known  a  priori  the  value  of  3  will  be  approximate.  As 

3  >  0  then  P(p)  >  0  for  all  p  and  is  a  continuous  function  of  p .  We 

must  now  determine  the  approximate  Legendre  coefficients,  This  is  done  by 

using  the  inversion  formula 
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2N+1 

2 


(36) 


0), 


1 

/  P(y)P  (y)dy. 
-1 


If  vi i>— 1  >  implying  usage  of  the  HG  phase  function,  none  of  the  to^  will 
exactly  match  the  exact  or  'true'  oo^.  Since  it  is  most  important  that  an 
approximate  phase  function  have  the  correct  single  scattering  albedo, 
w0,  we  redefine 


GO 

P(y)  =  P (y ) , 
GO 

0 


(37) 


and  therefore, 


go' 

N 


J  P(y)PN(y)dy 
-1 


GO 


GO 


or. 

N 


(38) 


This  guarantees  that  u>0  =  u0,  but  because  a,  g,  and  3  are  approximate, 
will  not  match  exactly.  As  long  as  the  HG  part  of  the  GHG  phase  is 

small,  the  mismatch  between  u>2  and  o>2  will  be  small. 

In  all  the  above  cases  it  is  necessary  to  calculate  only  two  of  the 

coefficients  in  the  Legendre  expansion  of  the  phase  function,  and 

do  or  to  and  to„.  This  may  be  done  directly  in  terms  of  the  a  and  b 
1  o  2  n  n 

coefficients  of  the  Legendre  series  occurring  in  Mie  theory  without  evaluating 
the  exact  phase  function  of  any  angle. 
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We  want  the  phase  function  incident  for  natural  light  to  satisfy 


/  P(y)df2  =  4ttcd  =  2t r  /  P(y)dy  =  4tt 
-1  -1  Cext 


(39) 


where  C  is  the  total  cross  section  for  scattering  and  extinction, 
respectively.  In  the  notation  of  van  de  Hulst4 


C 


scat 


-  y  /  F(0,<f>)df2. 
k 


(40) 


So,  if  we  put 


2tt 

P(y)=l?^  F(e,«d4»,  <41> 

o 


then 


/P(y)df2  =  K  k2  C 


scat ; 


(42) 


4H.  C.  van  de  Hulst,  1957,  Light  Scattering  by  Small  Particles,  John  Wiley  and 
Sons,  Inc.,  New  York  ~  ~~ 
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; 


or 


K  = 


(43) 


and  therefore 


P(y) 


k2  C 


ext 


2tt 

I  F(0,cf>)d<f>. 
o 


(44) 


Now  F(0,cj>)  = 


(0) sin 


<f>  + 


i2(0)cos  <p. 


hence 


P(y) 


k2  C 


ext 


[ijO) 


i2(0)]. 


\ 

2 

where  ilj2(0)  =  I sx  2(0)  I 


(45) 


and 


L-*» 

S  (0)  =  S  (0)  =  2  l  (2n+l)  (a  +b  )  . 

X  ^  **  „_-i  II  II 


(46) 
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Therefore, 


P(l) 


k2  C 


ext 


I  C2n+l)(a  +b  )|2  =  2 

n=l  n  n  k2  C 


ISjfO) 


(47) 


ext 


where 


C 


ext 


Re  [S1(0)] . 


(4  8) 


This  allows  direct  calculation  of  P(l)  in  terms  of  an  and  bn.  We  now  need 
expressions  for  and/or  io 2  in  terms  of  the  an  and  bn.  These  are  given  by 

Chu  and  Churchill. In  terms  of  the  Mie  coefficients  a  ,  b  ,  these 

n  n 

expressions  are: 


co. 


*Xxt  " 


^  t^?Slf)tVn*l+b„bn*l+an+lVb„.lV4  ^1)  <anVanV>- 


(49) 


5c.  M.  Chu  and  S.  W.  Churchill,  1955,  "Representation  of  the  Angular  Dis¬ 
tribution  of  Radiation  Scattered  by  a  Spherical  Particle,"  J  Opt  Soc  Am, 

45:958 

*There  appears  to  be  a  typographical  error  in  Chu  and  Churchill's  work  on  the 
bottom  of  page  961:  j+k-n  =  2r+l  should  read  j+k-n  =  2r-l. 
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4 


co 


2 


x2Q, 


ext 


00 


l 

n=l 


(n(n+l)-3)2C2n+l) 
n(n+l) (2n+3) (2n-l) 


*  * 

Re (a  a  +b  b  ) 
n  n  n  n 


15  n(n+3) 
2  (2n+3) 


"k  ie 

Re  (a  a  +b  b  ) 
n+2  n  n+2  n 


15  *  * 

+  Re(a  b  +b  .a  )] 
n+1  n+1  n  n+1  nJJ 


(50) 


where  x  =  2irr/X,  Qext  =  cext^irl‘  an<*  r  is  tlie  Particle  radius. 

RELATIONSHIPS  BETWEEN  SCATTERING  FRACTIONS  AND  PHASE  FUNCTIONS 


Among  the  objectives  of  the  work  done  here  was  the  conversion  of  the  single 
scattering  code  AGAUS  into  a  form  which  produced  the  types  of  data  required  by 
the  ACT  Battlefield/ Smoke  Obscuration  Model  while  retaining  a  variety  of 
options  previously  developed  for  AGAUS.  That  conversion  required  that  AGA.US 
be  given  the  additional  capabilities  of  predicting  (a)  attenuation 
coefficients  or  cross  sections  in  units  of  square  meters  per  milligram  of 
aerosol  material,  and  (b)  so-called  "scattering  fractions"  per  unit  aerosol 
mass.  In  order  to  avoid  possible  confusions  on  the  precise  relationships 
between  "scattering  fractions"  and  the  customary  quantities  used  in  Mie 
theory,  those  relationships  will  be  summarized  below  before  proceeding  to  a 
discussion  of  the  evaluation  of  AGAUSX. 

The  so-called  "phase  functions",  denoted  here  by  Pf(y)  and  "scattering 
fractions"  Sf(0)  are  both  derived  from  the  average  scattered  intensities 
(II  +  l2)/2  defined  in  Mie  theory.  The  two  quantities  are  related  to  one 

another  for  calculations  associated  with  a  single  wavelength  by  a  simple 
multiplicative  factor,  but  they  have  different  interpretations  and 
applications. 

Let  I(a,  m,  k,  9)  be  the  average  of  the  intensities  i^(0)  and  ^(Q)  for 
scattering  at  angle  0  from  a  sphere  whose  Mie-size  parameter  is  a  =  2irr/\,  and 
whose  complex  index  of  refraction  is  n  =  m-ik.  Furthermore,  let  n(r)dr  be  the 
relative  number  of  aerosol  particles  with  radii  between  r  and  (r+dr),  and  let 
Qgxt-Ca.mjk)  and  Qgca(a,m,k)  respectively  be  the  total  extinction  and 
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scattering  efficiency  factors  as  defined  by  van  de  Hulst.^ 
Now,  define  the  quantities  I,  Cexfc  and  Cgca  as  follows: 


C 


ext 


1  f  2  ,  ^ 

NT  J  nW 

T  r=o 


Qext (a>m>k)dr. 


(51) 


C 

sea 


1  r  2 

iC  J  1Tr 

T  r=o 


Qsca(a’m>k)dr> 


(52) 


1(6)  =  /  l(a,m,k,0)n(r)dr,  <53> 

T 

r=o 


where 


00 

NT  =  /  n(r)dr. 

r=o 


(54) 


^H.  C.  van  de  Hulst,  1957,  Light  Scattering  by  Small  Particles,  John  Wiley  and 
Sons,  Inc.,  New  York 
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The  "phase  function"  as  defined  in  program(s)  AGAUS  is  given  by 


Pf(6) 


2 

X 


TT  C 


ext 


1(9)  , 


(55) 


as  has  the  property  that 


2f  TT  C 

I  /  pf(0)d(j)  d (cos0)  =  4tt  =  4TT0J  ,  (56) 

cj)=o  0=0  r  0 

ext 


where  toQ  is  called  the  albedo  for  single  scattering" . 
The  "scattering  fractions"  Sf(0)  are  defined  by 


s£(0)  =  )  Nt  1(0), 


(57) 


and  S£(0)  has  the  property 


2Tr  TT 

/  I  S  (0)dcf)d(cos0)  =  C  *N 
(f)=o  0=o  sea  T 


(C  ) 

sea' total ; 


(58) 


wherein  (Csca)total  is  the  total  scattering  cross  section  per  unit  mass  of 
aerosol  material. 
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Thus,  it  will  be  seen  that 


s.(0)  =  (4-  c  j 

f  J  Mtt  T  ext^ 


P£(0) 


(59) 


The  ACT  Model  is  coded  with  the  assumption  that  (csca)total  as  derived  from 
Sj(6)  will  have  units  of  square  meters  per  milligram.  Programs  AGAUS9  and 
AGAUSX,  (see  below)  on  the  other  hand  are  coded  in  what  are  basically  CGS 
units.  Conversion  from  1(9)  to  S^(9)  in  the  ACT  normalization  therefore 

requires  a  unit-conversion  factor  for  from  pm  to  m,  and  for  from  cm  to 
-3  • 

m  ,  i.e.  , 


[A(Um)]2  NT(cm"3)  x  (10  6  j^)2 


x  106 


3 

cm 


3 

m 


(60) 


Thus,  a  factor  lO-6  is  required  to  convert  the  scattering  fractions  from  the 
internal  units  of  program  AGAUS  to  the  units  expected  by  ACT. 

VERIFICATION  OF  COMPUTER  CODES 

Evaluation  of  AGAUS9:*  Once  the  relationships  between  the  Mie  intensities,  ^ 
and  i  described  above  were  clearly  eluciated,  it  seemed  that  it  would  be  a 
straightforward  task  to  convert  AGAUS  from  the  computation  of  phase  functions 
to  the  calculation  of  scattering  fractions.  When  the  appropriate  conversions 
were  completed,  the  new  version,  AGAUS9,  was  run  using,  as  nearly  as  could  be 


*  AGAUS  9  is  the  parent  code  for  AGAUSX  and  they  basically  differ  only  in  the 
manner  in  which  the  number  of  radial  increments  for  integration  over  the  size 
distributions  are  chosen. 
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determined,  the  same  aerosol  model  for  which  results  were  found  in  the  ACT 
documentation.*  The  particular  model  used  was  a  log  normal  distribution  for 
white  phosphorus  (WP)  smoke  (r  =  0.37pm,  a  -  1.54pm),  200  particle  radii 
between  minimum  and  maximum  values  of  0.005  and  1.0pm,  respectively,  a  mass 
density  of  1.87  gm/cc,  a  particle  number  density  of  1.276345  x  10^  cm“^,  and 
complex  index  of  refraction  n  =  1.43  -  Oi.  Computations  were  performed  at  7 
wavelengths  (0.40,  0.45,  0.50,  0.55,  0.60,  0.65  and  0.‘70pm)  ,  and  the 
scattering  fractions  were  arithmetically  averaged  over  wavelength.  Some  of 
the  results  are  presented  in  table  1.  Only  a  cursory  examination  of  that 
table  is  needed  to  observe  that  the  AGAUS9  results  were  not  in  very  good 
agreement  with  those  found  in  the  ACT  documentation.  This  comparison  was 
unexpected,  because  the  basic  AGAUS  code  had  previously  been  found  to  give 
much  better  agreement  with  other  Mie  codes  than  is  seen  in  table  1. 

Because  thorough  reviews  of  the  program  listings  failed  to  shed  any  light  on 
the  reasons  for  such  large  apparent  discrepancies  it  was  decided  that 
comparisons  of  computations  made  at  a  single  wavelength  (rather  than  after 
averaging  over  several  wavelengths)  were  highly  desirable.  The  problem  there 
was  that  the  ACT  documentation  did  not  contain  any  single  wavelength 
results.  Since  that  avenue  for  testing  of  AGAUS 9  was  closed,  some  other 
approach  was  needed,  and  was  found  through  a  copy  of  the  MIE2  code. 

The  aerosol  model  used  to  generate  table  1  was  then  passed  through  both  AGAUS9 
and  the  MIE2  code  using  a  wavelength  of  0.40pm  only.  The  results  of  those  two 
runs  are  shown  in  table  2.** 

Examination  of  table  2  reveals  that  the  AGAUS 9  results  and  the  MIE2  results 
are  in  much  better  agreement  than  those  seen  in  table  1.  Whereas  table  2 
shows  discrepancies  even  in  the  first  digit  at  many  scattering  angles,  table  2 


*The  label  ACT  is  used  to  refer  to  the  document  entitled  "The  Effectiveness  of 
Obscuring  Smokes,"  by  Johnson,  Forney  and  Dolce,  and  unpublished  description 
of  the  JTOG/ME  smoke  obscuration  model  which  was  obtained  from  R.  B.  Gomez  of 
ASL. 


**The  results  of  MIE2  have  been  multiplied  by  a  factor  of  102  to  offset 
different  normalizations  used  in  ACT  and  AGAUS9. 
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TABLE  1 


COMPARISON  OF  SCATTERING  FRACTIONS  FROM  AGAUS9 
AND  ACT  REPORT  AS  AVERAGED  FOR  A  SEVEN  WAVELENGTH  RUN 


Scattering 

Angle 

Scattering  Fractions 

Scattering  Fractions 
ACT  Report 

0° 

0.5840074-02 

0.530822-02 

10° 

0.3385669-02 

0.294243-02 

20° 

0.1236144-02 

0.961893-03 

30° 

0.5525204-03 

0.419528-03 

O 

O 

0.2882904-03 

0.254840-03 

50° 

0. 1637546-03 

0.176892-03 

O 

O 

v£> 

0.9840574-04 

0.114931-03 

o 

O 

r-. 

0.6277324-04 

0.706178-04 

o 

o 

00 

0.4262488-04 

0.445657-04 

O 

o 

0.3097364-04 

0.309732-04 

100° 

0.2429897-04 

0.249208-04 

110° 

0.2086846-04 

0.237614-04 

120° 

0.2023009-04 

0.265251-04 

130° 

0.2267312-04 

0.330425-04 

140° 

0.3014750-04 

0.446135-04 

150° 

0.4592933-04 

0.637035-04 

160° 

0.6781825-04 

0.883152-04 

170° 

0.4512012-04 

0.678717-04 

180° 

0.6003613-04 

0.836141-04 
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TABLE  2 


COMPARISON  OF  RESULTS  FROM  RUNS  OF  AGAUS9 
AND  MIE2  AT  A  SINGLE  WAVELENGTH 


Scattering 

Angle 

Scattering  Fractions 

I*  (MIE2) 
avg  ' 

0° 

8.215439-03 

8.215442-03 

10° 

3.341243-03 

3.341238-03 

20° 

8.565565-04 

8.565530-04 

30° 

4.126651-04 

4.126640-04 

0 

O 

2.452253-04 

2.452247-04 

50° 

1.503114-04 

1.503111-04 

60° 

9.399925-05 

9.399906-05 

O 

O 

6.098411-05 

6.098398-05 

00 

o 

o 

4.167272-05 

4.167267-05 

O 

o 

3.019892-05 

3.019886-05 

100° 

2.333270-05 

2.333268-05 

110° 

1.957286-05 

1.957282-05 

120° 

1.854355-05 

1.854351-05 

130° 

2.096545-05 

2.096542-05 

140° 

2.942918-05 

2.942910-05 

150° 

4.947500-05 

4.947498-05 

160° 

8.643879-05 

8.643865-05 

170° 

6.208783-05 

6.208780-05 

180° 

7.644555-05 

7.644557-05 

2 

*The  MIE2  results  have  been  multiplied  by  10  to  normalize  them  in  the 
same  way  as  found  in  the  SOM  document. 
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shows  that  AGAUS9  and  MIE2  agree  to  at  least  give  and  often  six  significant 
digits — the  nominal  accuracy  which  is  used  in  terminating  Mie  series 
calculations  in  AGAUS9.  The  disagreements  found  in  table  1  should  be 
interpreted  as  differences  either  in  the  procedure  for  the  integration  over 
wavelength  or  the  fact  that  a  single  precision  version  of  MIE2  was  used  in  the 
generation  of  the  ACT  data. 

To  avoid  drawing  erroneous  conclusions  on  the  cross-agreements  between  results 
produced  by  AGAUS9  and  MIE2  on  the  basis  of  a  single  aerosol  model,  additional 
comparisons  were  made  using  a  quite  different  aerosol  model.  The  model  chosen 
for  this  "test"  was  a  "modified"  gamma  distribution: 


f(r)  =  Ara  exp(-BrY),  [for  MIE2] 


(61) 


with  A  =  1.0396  x  10®,  a  =  7.5,  B  =  333333  and  y  =  1*0.  The  input  data  used 
by  AGAUS9  are  the  particle  number  density  NT  and  mode  radius  r£,  rather  than  A 
and  B,  but  the  various  quantities  are  related  to  one  another  through: 


r 

c 


(— ) 

LyBJ 


1/y 


and 


N„ 


AB 


.(^1) 

Y 


r(— ) 

Y  J 


(62) 


For  this  comparison  the  minimum  and  maximum  particle  radii  were  taken  to  be 
0.01am  and  1.5pm,  respectively,  and  200  individual  values  of  particle  radii 
were  used.  The  run  was  made  at  a  wavelength  of  0.6pm  using  n  =  1.53  - 
0.006i.  Some  results  of  the  first  runs  of  this  type  will  be  found  in  table  3. 

Contrary  to  the  excellent  agreement  found  between  MIE2  and  AGAUS9,  the  data 
found  in  table  3  did  not  agree  as  well  as  had  been  expected  being  good  only  to 
about  4  digits.  Further  analysis  indicated  that  the  most  probable  source  of 
these  disagreements  lay  in  the  fact  that  AGAUS9  and  MIE2  did  not  choose  the 
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TABLE  3 


COMPARISON  OF  SCATTERING  FRACTIONS  CALCULATED  BY  MIE2 
AND  AGAUS9  FOR  A  MODIFIED  GAMMA  DISTRIBUTION 


Scattering 

Angle 

Scattering  Fractions 

I*  (MIE2) 
avg 

0° 

9.254421-08 

9.254768-08 

10° 

7.620680-08 

7.620913-08 

0 

O 

<M 

4.470319-08 

4.470374-08 

O 

O 

to 

2.156754-08 

2.156747-08 

40° 

1.038258-08 

1.038255-08 

50° 

5.626184-09 

5.626226-09 

60° 

3.415889-09 

3.415935-09 

O 

O 

r-. 

2.236736-09 

2.236776-09 

o 

o 

00 

1.555144-09 

1.555170-09 

O 

o 

1.149602-09 

1.149602-09 

100° 

9.113908-10 

9.114065-10 

110° 

7.876338-10 

7.876470-10 

120° 

7.656564-10 

7.656688-10 

130° 

8.657526-10 

8.657658-10 

140° 

1.074210-09 

1.074245-09 

150° 

1.193137-09 

1.193200-09 

160° 

1.134702-09 

1.134755-09 

170° 

1.502157-09 

1.502236-09 

180° 

1.969069-09 

1.969203-09 

*The  MIE2 

same  way 

2 

results  have  been  multiplied  by  10 

as  found  in  the  ACT  document. 

to  normalize  them 
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values  of  particle  radii  actually  used  in  the  calculation  in  the  same  way. 
Slight  alterations  to  remove  these  differences  were  then  made  in  AGAUS9,  and 
the  calculations  were  repeated,  yielding  the  data  presented  in  table  4. 

In  table  4  discrepancies  between  the  AGAUS9  and  MIE2  results  again  appear  only 
in  the  fifth  or  sixth  significant  digit. 

The  above  results  indicate  once  again  that  code  AGAUS9  is  quite  capable  of 
yielding  results  which  are  as  reliable  as  those  of  MIE2.  The  changes  in 
AGAUS9  nqeded  to  bring  about  the  level  of  agreement  seen  in  table  4  involved 
only  the  choice  of  the  values  of  the  200  radii  used,  and  illustrates  that  the 
method  of  choosing  the  radii  can  have  a  significant  effect  on  the  values  of 
the  scattering  fractions  even  when  a  relatively  large  number  of  radii  are 
used.  The  matter  will  be  discussed  further  in  a  subsequent  section  of  this 
document. 


EVALUATION  OF  PROGRAM  AGAUSX 

Introduction:  In  the  development  of  the  actual  coding  of  program  AGAUSX, 
continuous  testing  and  comparison  of  results  produced  by  the  new  codes  and 
those  yielded  by  other  single  scattering  codes  was  carried  out.  The  findings 
of  a  few  such  comparisons  will  be  given  next. 

The  supplementary  codes  used  for  evaluation  of  AGAUSX  were  the  codes  AGAUS9, 
and  the  code  MIE2.  In  order  to  avoid  the  problem  of  deciding  which  of  the 
three  codes  used  was  "correct"  in  the  event  that  no  two  agreed  well,  all  codes 
were  run  using  aerosol  models  for  which  independent  results  could  be  found  in 
Deirmendjian' s  book.^  The  specific  models  used  for  comparison  runs  were 
Deirmendjian' s  cloud  models  C.l  and  C.3  at  a  wavelength  of  0.7pm.  Additional 
comparison  runs  were  made  at  10.6pm  using  a  simulated  hygroscopic  smoke  model^ 


^D.  Deirmendjian,  1969,  Electromagnetic  Scattering  on  Spherical 
Polydispersions,  American  Elsevier  Publishing  Co.,  Inc.,  New  York 

7 

A.  Miller,  R.  C.  Shirkey,  E.  Gemoets  and  G.  H.  Goedecke,  1978,  Investigations 
on  the  Prediction  of  Infrared  Transmission  and  Emission  by  Clear  and  Aerosol 
Laden  Atmospheres"  NMSU  Department  of  Physics ,  Final  Report  for  Contract 
DAEA18-77-C-00003 
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TABLE  4 


COMPARISON  OF  SCATTERING  FRACTIONS  PRODUCED  BY  MIE2  AND  A  SPECIAL  VERSION* 
OF  AGAUS9  FOR  A  MODIFIED  GAMMA  DISTRIBUTION 


Scattering 

Angle 

Scattering  Fractions 

I*  (MIE2) 
avg'' 

0° 

9.254772-08 

9.254768-08 

1— » » 
O 
o 

7.620918-08 

7.620913-08 

N) 

O 

o 

4.470380-08 

4.470374-08 

04 

o 

o 

2.156751-08 

2.156747-08 

o 

O 

1.038257-08 

1.038255-08 

o 

O 

LO 

5.626233-09 

5.626226-09 

60° 

3.415938-09 

3.415935-09 

O 

O 

2.336777-09 

2.336776-09 

o 

O 

00 

1.555171-09 

1.555170-09 

o 

o 

1.149623-09 

1.149622-09 

100° 

9.114071-10 

9.114065-10 

110° 

7.876480-10 

7.876470-10 

120° 

7.656697-10 

7.656688-10 

130° 

8.657668-10 

8.657658-10 

140° 

1.074245-09 

1.074245-09 

150° 

1.193200-09 

1.193200-09 

160° 

1.134755-09 

1.134755-09 

170° 

1.502236-09 

1.502236-09 

180° 

1.969203-09 

1.969203-09 

Particle  radii  selected  the  same  way  as  in  MIE2. 

*  2 
The  MIE2  results  have  been  multiplied  by  10  to  normalize  them  in  the 

same  way  as  found  in  the  SOM  document. 
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and  several  values  of  relative  humidity.  The  latter  comparisons  involved  only 
program  AGAUS9  and  AGAUSX.  The  objectives  of  all  comparisons  were  not  only  to 
show  that  AGAUSX  was  working  properly,  but  also  to  determine  whether  or  not  it 
would,  in  fact,  offer  significant  reductions  in  overall  running  time. 

Water  Cloud  Model  Comparisons:  Comparison  runs  of  codes  AGAUS9,  AGAUSX  and 
MIE2  for  cloud  models  C.l  and  C.3  were  made  using  the  same  range  of  radii 
quoted  for  Deirmend jian* s  tables  T.  36  and  T.30.  In  order  to  assure  that  any 
conclusions  drawn  about  relative  computation  times  would  be  valid  in  the 
context  of  ACT  usage,  scattering  fractions  were  calculated  at  2°  increments 
between  scattering  angles  of  0°  and  180°  in  these  runs.  AGAUSX  was  also  run 
for  model  C.l  using  several  different  values  of  the  convergence  level  DELTA. 

Table  5  presents  some  of  the  results  obtained  for  cloud  model  C.l.  The 
scattering  fraction  entries  for  the  "Deirmend jian"  column  were  obtained  by 
averaging  his  Pj  and  P2,  and  multiplying  by  the  appropriate  conversion 
factor.  The  AGAUSX  run  listed  used  DELTA  =  0.01.  Examination  of  the 
tabulated  data  shows  that  none  of  the  three  runs  yields  exactly  the  same 
extinction  coefficient  that  was  printed  by  Deirmend jian,  but  all  were  within 
0.3  percent  of  his  results.  Unfortunately,  Deirmend jian' s  table  T.36  did  not 
include  angles  at  2°  increments,  so  only  forward  and  backward  scattering 
fractions  could  be  inferred  accurately.  It  will  be  seen  that  all  three  codes 
gave  scattering  fractions  which  agreed  well  in  the  0=0°  direction,  but  the 
results  are  obviously  not  identical.  The  variations  are,  however,  quite  a  bit 
larger  at  other  scattering  angles,  although  the  variations  are  not  terribly 
large  when  their  magnitudes  are  compared  to  that  of  the  forward  scattering. 
For  this  particular  model,  it  will  be  noted  that  AGAUSX  convergence  criterion 
led  to  its  use  of  more  radii  than  were  used  with  the  other  codes.  It  will 
also  be  seen  that  MIE2  handled  the  40  radii  faster  than  did  AGAUS9.  MIE2's 
speed  advantage  over  AGAUS9  probably  results  from  the  fact  that  MIE2  used  non- 
uniform  spacing  of  particle  radii — a  procedure  which  leads  to  a  significantly 
smaller  number  of  Mie  calculations  which  must  be  done  ar  large  Mie  size- 
parameters. 
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COMPARISONS  OF  RESULTS  AND  CPU  TIMES  FOR  THREE  SINGLE  SCATTERING  CODES 
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*N0TE :  AGAUS9  and  AGAUSX  scattering  fractions  at  0°,  when  converted  to  phase  functions  have  values 
of  1683,  as  compares  to  Deirmendjian' s  values  of  1680. 


It  should  be  noted  that  the  version  of  AGAUSX  used  in  generating  results  shown 
in  table  5  performed  its  "convergence  test"  on  the  total  aerosol  volume  and 
not  on  the  scattering  fractions.  (Results  obtained  by  checking  on  both  volume 
and  the  scattered  intensity  at  90°  will  be  given  below.) 

Although  AGAUSX  appears  to  offer  no  computation  speed  advantage  if  the  user 
knows  how  many  radii  are  "enough",  it  does  have  the  useful  feature  (not  found 
either  MIE2  or  AGAUS9)  of  giving  the  user  some  idea  of  how  many  radii  were 
really  required  to  achieve  some  minimum  level  of  confidence  in  its  results. 

By  repeating  runs  of  AGAUSX  and  model  C.l  using  various  values  of  DELTA,  it 
has  been  possible  to  learn  more  about  the  sensitivities  of  the  extinction 
coefficient  and  scattering  fractions  to  the  number  of  radii  used  in  a 
calculation.  Table  6  contains  the  results  of  a  few  runs  of  that  type.  The 
first  four  entries  represent  runs  in  which  the  quantity  which  was  tested  for 
convergence  was  the  total  volume  of  the  aerosol  particles.*  The  fifth  and 
sixth  entries  are  the  results  of  runs  in  which  testing  was  done  on  both  volume 
and  the  scattered  intensity  at  90°;  note  that  the  sixth  entry  represents  a  run 
at  45°  increments  rather  than  2°  increments. 

Examination  of  table  6  shows  that  the  value  of  the  extinction  coefficient  is 
not  nearly  as  sensitive  to  the  choice  of  particle  radius  values  as  are  the 
scattering  fractions  at  fairly  large  scattering  angles,  especially  near  90°. 
The  90°  results  show  "oscillations"  of  as  much  as  15  percent  as  the  number  of 
radii  changes  through  the  set  {26,  50,  66,  98,  130,  and  513}. 

It  should  also  be  noted  that  AGAUS9  and  MIE2  runs  shown  in  table  5  used  only 
40  radii,  while  Deirmendjian  evidently  used  440  radii.  Even  with  that  vastly 
larger  number  of  radii,  Deirmendjian' s  results  are  not  so  different  from  the 
present  results  to  warrant,  in  many  applications,  an  additional  increase  in 
computation  time  of  a  factor  of  eleven. 


*Volume  was  used  rather  than  extinction  cross-section  because  its 
dependence  appeared  to  make  it  a  "more  sensitive"  test  quantity. 
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COMPARISON  OF  AGAUSX  RESULTS  OBTAINED  WITH  FOUR  CONVERGENCE  TEST  LEVELS 
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A  similar  set  of  comparative  results  obtained  for  cloud  model  C.3  are 
presented  in  table  7.  For  this  model,  AGAUS9  and  MIE2  were  run  with  the  same 
number  of  radii  as  was  used  by  Deirmend jian.  This  example  illustrates  the  way 
in  which  AGAUSX  can  offer  definite  decreases  in  CPU  time  over  codes  requiring 
the  user  to  make  an  a  priori  decision  as  to  the  number  of  radii  needed.  Using 
only  66  radii,  AGAUSX  produced  an  extinction  coefficient  to  within  better  than 
1  percent  of  Deirmend jian' s  value  (with  which  AGAUS9  and  MIE2  results  are 
identical  to  the  four  digits  given  by  Deirmend jian) .  AGAUSX,  (with  DELTA  = 
0.01),  however,  ran  nearly  four  times  faster  than  either  AGAUS9  or  MIE2. 
Furthermore,  the  AGAUSX  run  at  DELTA  =  0.01  yielded  scattering  fractions 
which,  at  worst,  differed  from  those  found  by  AGAUS9  and  MIE2  by  the  order  of 
10  percent. 

Cloud  models  such  as  C.l  and  C.3  used  at  short  wavelengths  ( s*  visible)  are 
severe  tests  of  single  scattering  codes  because  they  typically  involve  some 
rather  large  Mie  size-parameters  (a  110  for  the  C.l  model).  Calculations 
performed  at  infrared  wavelengths  for  smaller  particles  than  found  in  typical 
clouds  are  not  quite  so  demanding.  That  fact  is  illustrated  partially,  by 
tables  8  and  9,  which  were  compiled  using  a  hypothetical  model  for  a 
hygroscopic  smoke  having  the  particle  size-distribution  of  white  phosphorous 
smoke.  The  model  used  here  was  denoted  as  model  A'  in  an  earlier  document.7 
Table  8  shows  extinction  coefficients  obtained  from  AGAUS9  (using  100  radii), 
and  AGAUSX  (using  DELTA  =  0.01)  at  seven  different  values  of  relative 
humidity.  It  also  presents  the  results  of  a  run  of  a  special  version  of 
AGAUSX  (labeled  AGAUSXL)  in  which  the  Mie  routine  D0WN42  was  used  instead  of 
routine  MIEGX.  It  will  be  seen  that  AGAUSX  (and  AGAUSXL)  needed  only  34  radii 
to  reach  extinction  coefficients  within  0.2  percent  of  those  found  by  AGAUS9 , 
but  the  AGAUS9  run  used  340  percent  more  computer  time  than  did  AGAUSX. 
AGAUSXL' s  results  were  quite  close  to  those  of  AGAUSX,  but  the  relative 
slowness  of  DCWN42  brought  the  total  time  to  treat  34  radii  up  to  nearly  the 
same  time  needed  by  AGAUS9  to  treat  100  radii. 


A.  Miller,  R.  C.  Shirkey,  E.  Gemoets  and  G.  H.  Goedecke,  1978,  Investigations 
on  the  Prediction  of  Infrared  Transmission  and  Emission  by  Clear  and  Aerosol 
Laden  Atmospheres,  NMSU  Department  of  Physics,  Final  Report  for  Contract 
DAEA18-77-C-00003 
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COMPARISONS  OF  RESULTS  AND  CPU  TIMES  FOR  THREE  SINGLE  SCATTERING  CODES 
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(b)  The  quantity  tested  for  convergence  was  the  total  aerosol  volume. 


TABLE  8 


COMPARISONS  OF  RESULTS  FROM  PROGRAMS  AGAUS9  AND  AGAUSX 

FOR  SMOKE  MODEL  A'  AT  A  =  10. 6pm  (Relative  humidity  =  0%) 

Extinction  Coefficients  (per  km) 

AGAUSX  AGAUSXL 

Relative  Humidity  AGAUS9  (34  radii)  (34  radii) 

(percent) _  (100  radii)  _ [Delta  -  0.01) _  %Difference * 


0 

0. 1838 

0.1835 

0.1835 

0.16 

75 

0.2801 

0.2797 

0.2797 

0.18 

80 

0.3020 

0.3016 

0.3016 

0.17 

85 

0.3340 

0.3335 

0.3335 

0.27 

90 

0.3887 

0.3882 

0.3882 

0.15 

95 

0.5350 

0.5343 

0.5343 

0.15 

— 

99 

1.9295 

1.9268 

1.9268 

0.15 

TOTAL 

CPU  TIME 

34.34  sec 

9.83  sec 

34.98  sec 

71.4 

*NOTE : 

Percent 

differences  are: 

| (AGAUS9  - 

AGAUSX) /AGAUS9|  x  100 
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TABLE  9 


COMPARISON  OF  AGAUSX  RESULTS  FOR  FOUR  DIFFERENT  CONVERGENCE  LEVELS 
AND  A  SIMULATED  SMOKE  MODEL  -  A  =  10.6ym  (Relative  humidity  =  0%) 


Convergence  Level 
(DELTA) 

No.  of  Radii  Used 

Kext<km 

.05 

18 

1.8277 

.01 

34 

1.8353 

.005 

50 

1.8345 

.001 

130 

1.8377 

Convergence  Level 
(DELTA) 

O 

O 

11 

CD 

Scattering  Fractions 

0  =  90° 

6  =  180' 

.05 

1 . 5574E-6 

6. 9366E-7 

1 . 2326E-6 

.01 

1 . 5610E-6 

6. 9555E-7 

1 . 2365E-6 

.005 

1 . 5605E-6 

6. 9533E-7 

1 . 2361E-6 

.001 

1 . 5621E-6 

6 . 9613E-7 

1 . 2376E-6 
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Table  9  demonstrates  how  the  user's  choice  of  DELTA  affected  runs  of  the 
simulated  smoke  model  at  a  relative  humidity  of  0  percent.  For  that  model, 
changes  in  the  number  of  radii  from  18  to  130  (a  factor  of  seven)  did  not 
change  either  the  extinction  coefficient  or  the  scattering  fractions  by  more 
than  0.3  percent. 

Discussion:  One  of  the  decisions  which  had  to  be  made  for  program  AGAUSX  was 
that  of  determining  just  which  "quantity"  or  set  of  quantities  should  be  used 
in  the  test  for  convergence  to  within  the  level  DELTA.  The  one  which  was 
finally  chosen  was  the  total  volume  of  the  aerosol  particles.  The  volume  was 
chosen  because,  being  dependent  on  the  cube  of  the  radii,  it  might  be  a  little 
more  conservative  choice  than  the  extinction  coefficient,  although,  as  seen  in 
tables  5  and  6,  it  is  a  less  conservative  test  quantity  than  the  scattered 
intensity  at  9  =  90°.  Present  experience  is  limited,  but  the  use  of  the 
volume  seems  to  be  a  good  compromise  unless  highly  precise  scattering 
fractions  are  needed.  If  the  latter  situation  arises  in  a  particular 
application,  users  of  the  code  should  have  no  difficulty  in  coding  in  their 
own  preferred  test  quantity.  The  quantities  which  are  available  at  present 
are  extinction,  scattering,  backscattering  cross  sections,  and  the  average 
intensities  at  various  angles.  Careful  study  of  the  source  listing  for 
subroutine  AGXPT2  will  reveal  how  such  changes  can  be  made. 

A  major  surprise  which  was  revealed  by  these  studies  was  the  fact  that  the 
code  MIE2  used  less  computer  time  to  handle  cloud  models  C.l  and  C.3  than 
AGAUS9  needed  (for  the  same  number  of  radii).  The  reason  for  the  surprise  was 
that  the  basic  Mie  routine  DCWN42  used  by  MIE2  was  known  to  be  usually 
appreciably  slower  than  routine  MIEGX.  It  was  in  fact,  that  difference  in 
speed  which  caused  AGAUSXL  (see  above,  and  table  8)  to  require  more  than  three 
times  as  long  to  treat  the  smoke  model  than  AGAUSX  needed.  The  reason  for 
MIE2's  evident  computation  time  advantage  over  AGAUS9  has  been  alluded  to 
above,  and  appears  to  be  the  result  of  the  differences  in  the  ways  in  which 
the  values  of  radii  (at  which  the  Mie  calculations  are  to  be  done)  are  chosen 
in  the  two  codes.  In  MIE2,  the  interval  between  successive  values  of  particle 
radius  is  doubled  as  each  additional  value  is  assigned  (at  least  for  the 
generalized  Khirgian-Mazin  distribution),  while  AGAUS9  uses  a  constant 
increment.  The  MIE2  method  therefore  uses  far  fewer  large  values  of  radius 
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than  the  AGAUS9  method  does.  That  difference,  and  the  fact  that  the  time 
required  within  the  Mie  subroutines  D0WN42  and  MIEGX  increases  rapidly  with 
the  size  to  wavelength  ratio,  appears  to  account  for  MIE2's  (small)  running 
time  advantage  over  AGAUS9.  That  difference  in  how  the  radii  are  chosen 
definitely  accounts  for  the  disagreements  between  the  MIE2  and  AGAUS9  values 
for  the  scattering  fractions  for  cloud  models  C.l  and  C.3  (that  disagreement 
is  quite  appreciable  in  table  5  in  some  cases).  The  latter  conclusion  has 
been  verified  by  direct  runs  with  a  version  of  AGAUS9  which  uses  the  MIE2 
method  for  choosing  the  radii  (the  AGAUS9/MIE2  disagreements  moved  out  to  the 
fifth  or  sixth  digit).  The  running  time  on  a  UNIVAC  1108  of  the  special 
version  of  AGAUS9  was  240  seconds  as  compared  to  266  seconds  for  MIE2, 
demonstrating  that  the  use  of  fewer  large  radii  significantly  decreases 
running  time. 

The  questions  raised  above  are,  however,  made  somewhat  academic  by  AGAUSX  for 
many  applications  requiring  single  scattering  calculations  because  of  AGAUSX' s 
clear  speed  advantage  in  handling  distributions  for  which  no  a  priori  guide  as 
to  how  many  radii  will  be  "enough"  is  available. 

SOME  EFFECTS  OF  THE  NIMER  OF  SIZE  INTERVALS  USED  IN  AGAUSX 

In  the  early  development  of  AGAUSX,  it  was  tentatively  assumed  that 
computation  times  might  be  substantially  reduced  by  breaking  the  total 
ranges  of  radius  values  to  be  treated  into  more  than  one  "size-interval"  and 
then  performing  the  halving  calculations  and  tests  separately  within  each  size 
interval.  That  assumption  was  based  upon  the  belief  that  the  use  of  a  single 
interval  could  result  in  unneeded  computations  in  some  size  ranges  brought 
about  by  slow  convergence  in  other  ranges.  Consequently,  all  types  of 
distribution  functions  which  were  known  to  be  "peaked"  were  split  into  two 
size  intervals:  (a)  one  with  radii  smaller  than  the  one  at  which  the 

distribution  was  a  maximum,  and  (b)  one  containing  all  larger  values  of 
particle  radii.  In  some  of  the  test  runs,  however,  it  was  found  that  the 
interval  containing  the  smaller  radii  required  more  halvings  than  the  other 
interval,  but  contributed  only  a  small  fraction  of  the  total  extinction 
effects — thereby  wasting  computation  time. 
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The  effect  of  using  just  a  single  size  interval,  instead  of  two  intervals,  has 
been  briefly  explored  using  a  special  version  of  AGAUSX,  and  Deirmend jian' s 
cloud  models  C.l  and  C.3  at  X  =  0.7am.  Computations  were  made  with  several 
different  values  for  the  convergence  level  (A).  Table  10  compares  some  of  the 
results  found  with  the  special  version  (AGAUS10S)  and  the  regular  version. 
The  scattering  fraction  computations  were  made  at  2°  increments  in  most  cases 
shown  in  the  table,  and  the  column  labeled  "NI"  explicitly  shows  the  number  of 
size  intervals  used.  Also  shown  and  labeled  "REF  CASE",  are  the  results 
obtained  with  an  AGAUSX  run  which  used  all  513  available  values  of  particle 
radii. 

The  results  seen  in  table  10  indicate  that  the  special  version  (AGAUS10S) 
required  only  one-third  to  one-half  as  much  computation  time  to  satisfy  the 
different  convergence  levels  as  AGAUSX  needed.  It  thus  appears  that  the 
initial  assumption  that  two  intervals  would  be  preferable  to  one  interval  may 
be  false.  However,  a  closer  examination  of  the  table  shows  that  AGAUS10S 
yielded  substantially  large  scattering  fractions  (away  from  0°)  than  did 
AGAUSX,  with  the  latter  code's  values  being  closer  to  the  "REF  CASE"  values  at 
all  choices  of  A. 

About  the  only  conclusions  which  can  be  drawn  at  this  point  are  that  the  NI  = 
1  version  is  by  far  the  most  efficient  when  only  kgxt  is  of  interest,  and  that 
the  NI  =  2  version  is  preferable,  in  spite  of  its  slower  speed,  if  high 
accuracies  are  needed  for  the  scattering  fractions  at  large  scattering 
angles.  From  another  vantage  point,  however,  the  discrepancies  between  the  NI 
=  1  and  NI  =  2  scattering  fractions  at  0  =  90°  and  9  =  180°  are  an 
insignificant  fraction  of  the  0=0  values. 

A  few  similar  comparisons  for  Deirmend jian' s  cloud  model  C.3  at  X  =  0.7pm  are 
given  in  table  11. 

Comparison  of  Results  for  the  Analytic  Phase  Functions:  Ttoo  models  have  been 
chosen  to  compare  the  two  analytic  phase  functions,  Henyey-Greenstein  and 
Goedecke-Henyey-Greenstein  with  the  'true'  or  computed  phase  function.  These 
models  and  selected  wavelengths  at  which  comparisons  were  made  are:  simulated 
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TABLE  10 


COMPARISON  OF  RESULTS  PRODUCED  FOR  CLOUD  MODEL  C.l  (A  =  0.7ym) 

BY  SEVERAL  VARIANTS  OF  PROGRAM  AGAUSX  AND  DIFFERENT  CONVERGENCE  LEVELS  (A) 


_ Scattering  Fractions 

DELTA  *  NI  NRADI  K  „  CPU  O73  9lT  ISTT 

ext 


a 

2 

26 

17.47 

45.28 

2.4521 

5 . 2019E-5 

7.4914E-4 

0.05 

b 

1 

9 

16.70 

21.06 

2.2353 

6 . 9736E-5 

1.29463-3 

c 

1 

9 

16.70 

note-d 

2.2353 

6. 9736E-5 

1 . 2946E-3 

a 

2 

50 

16.72 

82.81 

2.2404 

4 . 8011E-5 

8 . 3856E-4 

0.01 

b 

1 

17 

16.77 

38.72 

2.2438 

7 . 0311E-5 

1 . 0222E-3 

c 

1 

129 

16.74 

note-d 

2.2386 

5 . 6969E-5 

9. 1422E-3 

0.005 

a 

2 

66 

16.74 

127.93 

2.2395 

5 . 6570E-5 

8 . 7430E-4 

b 

1 

17 

16.77 

37.78 

2.2438 

7. 0311E-5 

1 . 0222E-3 

130 

16.73 

245.13 

2.2370 

5 . 4473E-5 

8 . 2393E-4 

33 

16.76 

75.71 

2.2393 

6 . 8529E-5 

9. 3237E-4 

"Reference  c 

Case" 

2 

513 

16.75  note-d 

2.2411 

5 . 5880E-5 

8 . 8647E-4 

Deirmendj ian 

1 

440 

16.73 

2.2370 

8.457E-4 

*NOTES: 


(a)  "Normal"  AGAUSX  convergence  tests  on  aerosol  volume;  2  size  intervals. 

(b)  AGAUS10S;  convergence  tests  on  extinction  coefficient;  1  size 
interval . 

(c)  AGAUS10S;  convergence  tests  on  extinction  coefficient  and  scattered 
intensity  at  0  =  90°;  number  of  size  intervals  given  in  column 
labeled  "NI". 

(d)  45°  angular  increment  runs;  CPU  times  not  directly  comparable  to 
other  runs  made  at  2°  increments. 
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TABLE  11 


COMPARISON  OF  AGAUSX  RESULTS  FOR  CLOUD  MODEL  C.3  (A  =  0.7ym)  USING  ONE 
AND  TWO  SIZE  INTERVALS,  SEVERAL  CONVERGENCE  LEVELS  AND  CONVERGENCE  TESTS 
ON  EXTINCTION  COEFFICIENT  AND  SCATTERING  FRACTIONS  AT  90° 


Scattering  Fractions 


DELTA 

NI 

NRADI 

K  * 
ext 

cpu(a;) 

0° 

90° 

180° 

0.1 

-- 

-- 

-- 

-- 

— 

_ 

_  _ 

1 

33 

3.042 

34.81 

5.4626E-2 

2 . 2433E-5 

2 . 5350E-4 

0.05 

-- 

-- 

-- 

-- 

— 

_ 

—  — 

1 

65 

68.71 

5.4583E-2 

2 . 2435E-5 

2 . 3498E-4 

2 

66 

71.49 

5 . 3569E-5 

2 . 1668E-5 

1 . 7792E-4 

1 

65 

3.041 

72.67 

5 . 4583D-2 

2 . 2435E-5 

2 . 3498E-4 

0.005 

-- 

-- 

-- 

— 

_ 

-  _ 

1 

65 

3.041 

72.06 

5 . 4583E-2 

2 . 2435E-5 

2 . 3498E-4 

0  001 

2 

130 

3.015 

5 . 3796E-2 

2 . 2180E-5 

1 . 9798E-4 

1 

129 

3.041 

16. 86b 

5 . 4585E-2 

2 . 1382E-5 

2 . 3500E-4 

_ k _ 

.0001 

1 

513 

62.83“ 

5 . 3985E-2 

Deirmendj ian 

280 

3.021 

5.396E-2 

_ 

2.025E-4 

NOTES:  (a)  CPU  times  in  seconds,  using  2°  angular  increments. 

(b)  Not  comparable  to  other  runs  due  to  use  of  45°  increments. 

(c)  Cases  for  which  dashes  ( - )  appear  were  not  run  for  NI=2 

and  convergence  tests  on  both  K  and  the  scattered  intensity 
at  0  =  90°. 
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white  phosphorus^  at  1.06  and  10.6  microns;  fog  model  WSMRF2  at  wavelengths  of 
2.5,  10.0  and  11.0  microns;  the  parameters  for  each  of  the  two  fog  models  are 
tabulated  below. 


WSMRF2 


Radius,  Minimum 

0. 4pm 

Radius,  Maximum 

1 2.  0pm 

Radius,  Mode 

4.  0pm 

Alpha 

6.0 

Gamma 

1.0 

Particle  Density 

100. 0cm' 

Values  of  the  analytical  phase  functions  HG  and  GHG,  along  with  the  computed 
phase  functions  are  presented  graphically  in  figures  1  through  5.  A  visual 
inspection  of  these  figures  shows  that  the  GHG  analytic  phase  function  matches 
the  computed  phase  function  far  better  than  the  HG  phase  function  near  zero 
degrees,  and  that  the  GHG  function  apparently  provides  a  better  overall  fit  to 
the  computed  phase  function  than  the  HG  function.  Because  the  scales  of  the 
graphs  usually  do  not  permit  sufficient  resolution  near  180° ,  the  values  of 
the  computed  and  analytic  phase  functions  are  presented  in  table  12  at 
selected  angles.  Table  12  also  presents  the  differences  between  the  computed 
and  analytic  phase  functions  at  the  selected  angles,  and  also  the  root  mean 
square  error  for  the  points  considered. 

Discussion:  Inspection  of  table  12  shows  that  the  GHG  analytic  phase  function 
is  superior  to  the  HG  analytic  phase  function  at  angles  near  zero  degrees. 


7 A.  Miller,  R.  C.  Shirkey,  E.  Gemoets  and  G.  H.  Goedecke,  1978,  Investigations 
on  the  Prediction  of  Infrared  Transmission  and  Emission  by  Clear  and  Aerosol 
Laden  Atmospheres,  NMSU  Department  of  Physcis,  Final  Report  for  Contract 
DAE Al 8-77 -C-00003 
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PHASE  FUNCTION 


0  45  90  135  180 

ANCLE 

(DEGREES) 


Figure  1.  Comparison  of  analytic  phase  functions  HG  and  GHG 
with  computed  phase  function  for  simulated  white 
phosphorous  at  1.06ym. 


ANGLE 

(DEGREES) 


Figure  2.  Comparison  of  analytic  phase  functions  HG  and  GHG  with 
computed  phase  function  for  simulated  white  phosphorous 
at  10.6ym. 
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PHASE  FUNCTION  PHJ!!?.E..f_l)NC.T_l0N 


Figure  4 


Comparison  of  analytic  phase  functions  HG  and  GHG 
with  computed  phase  functions  for  fog  model  WSMEF2 
at  2.5pm. 


Comparison  of  analytic  phase  functions  HG  and  GHG 
with  computed  phase  functions  for  fog  model  WSMRF2 
at  10.0pm. 


PHRSE  FUNCTION 

(ARBITRARY  UNITS) 


Figure  5.  Comparison  of  analytic  phase  functions  HG  and  GHG  with 
computed  phase  function  for  fog  model  WSMRF2  at  ll.Oum. 


TABLE  12 


COMPARISON  OF  COMPUTED  AND  ANALYTIC  PHASE  FUNCTION  VALUES  AT  SELECTED  ANGLES  FOR 
FOG  MODEL  WSMRF2  AND  SIMULATED  WHITE  PHOSPHOROUS  AT  VARIOUS  WAVELENGTHS 


Wavelength 
2. 5iim 


Model  WSMR  F2 
Phase  Function  Values* 


Angle(°) 

True  t 

HG 

HG-GG 

1.51 

134.80 

68.17 

133.20 

66.63 

1.60 

25.08 

3.59 

3.56 

2.56 

.03 

1.30 

50.66 

.49 

.55 

.43 

-.06 

.06 

74.26 

.13 

.20 

.18 

-.07 

-.05 

99.84 

.06 

.10 

.11 

-.04 

-.05 

125.41 

.09 

.06 

.08 

.03 

.02 

150.98 

.09 

.05 

.07 

.04 

.02 

174.58 

.17 

.05 

.07 

.12 

.10 

178.49 

.23 

.05 

.07 

.18 

.16 

22.21 

.64 

Wavelength 

Model 

WSMR  F2 

1 0 . Oym 

Phase  Function  Values* 

1.51 

14.12 

60.63 

12.86 

-46.51 

1.26 

25.08 

4.32 

2.07 

3.97 

2.25 

.35 

50.66 

.38 

.30 

.68 

.08 

-.30 

74.26 

.08 

.11 

.09 

-.03 

-.01 

99.84 

.03 

.05 

.00 

-.02 

.03 

125.41 

.02 

.04 

.00 

-.02 

.02 

150.98 

.02 

.03 

.00 

-.01 

.02 

174.58 

.02 

.02 

.00 

.00 

.02 

178.49 

.02 

.02 

.00 

.00 

.02 

15.52 


.45 
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However,  as  one  approaches  180°  the  HG  function  sometimes  provides  a  better 
fit  to  the  computed  phase  function.  Since  for  most  of  the  distributions 
examined  the  backscatter  is  far  less  than  the  forward  scattering,  it  becomes 
important  to  analytically  fit  the  computed  phase  function  more  closely  near 
the  forward  direction  of  scattering  relative  to  the  backscatter  direction. 
Additionally  since  the  HG  analytic  function  fits  the  computed  phase  function 
best  near  180°  and  the  GHG  analytic  function  is  partially  comprised  of  the  HG 
function  near  180°,  we  again  find  that  the  GHG  analytic  function  will  fit  well 
near  180°.  The  root  mean  square  error  is  also  presented  for  the  difference 
between  the  computed  and  analytic  phase  functions,  and  again  indicates  that 
the  GHG  analytic  phase  function  provides  a  better  fit.  This  is  particularly 
so  if  the  particular  regime  under  consideration  is  near  the  incident  (forward) 
direction. 


STUDIES  OF  IR  EXTINCTION  IN  HYGROSCOPIC  AEROSOLS 

Introduction:  In  a  previous  document^  some  results  of  the  modeling  of  IR 

extinction  in  hygroscopic  aerosols  were  reported.  Since  then  additional 
computations  have  been  made  to  extend  the  range  of  relative  humidities  beyond 
that  used  previously  and  to  examine  the  effects  of  "hysteresis"  in  the 
condensation  and  evaporation  of  droplets  formed  upon  hygroscopic  condensation 
nuclei.  The  new  calculations  were  carried  out  using  the  simulated  white 
phosphorous  smoke  model  A'  defined  in  the  report  referenced  above. 

Cbmputations  were  carried  out  with  an  early  version  of  program  AGAUS9  at 
wavelengths  of  0.55,  1.06,  and  3.6um.  Twelve  values  of  relative  humidity  were 
used,  and  mass  accretion  coefficients  were  taken  from  Hanel^  for  conditions  of 


^ A.  Miller,  R.  C.  Shirkey,  E.  Gemoets  and  G.  H.  Goedecke,  1978,  Investigations 
on  the  Prediction  of  Infrared  Transmission  and  Emission  by  dear  and  Aerosol 
Laden  Atmospheres,  NMSU  Department  of  Physics,  Final  Report  for  Cbntract 
DAEA18-77-C-00003 

®G.  Hanel,  1976,  "The  Properties  of  Atmosphere  Particles  as  Functions  of  the 
Relative  Humidity  at  Thermodynamic  Equilibrium  with  the  Surrounding  ffoist 
Air,”  Ad van  Geophys ,  Editors  H.  E.  Landsberg  and  J.  Van  MLeghem,  Academic 
Press,  New  York 
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Results :  The  data  generated  by  the  new  computations  are  presented  in 

graphical  form  in  figures  6  through  8.  Figure  6  shows  the  absolute  extinction 
coefficients  as  a  function  of  relative  humidity  at  X  -  0.55pm  and 
X  =  1.06pm.  The  arrowheads  indicate  the  direction  in  which  relative  humidity 
is  changing;  an  arrowhead  tilted  toward  the  right  represents  increasing 
relative  humidity,  and  an  arrowhead  tilted  toward  the  left  represents 
decreasing  relative  humidity.  Figure  6  clearly  shows  the  "hysteresis”  effect 
mentioned  above.  It  will  be  seen  that  for  relative  humidities  between  about 
60  and  75  percent,  the  value  of  the  extinction  coefficient  may  vary  by  a 
factor  of  two  between  conditions  of  rising  and  falling  relative  humidity.  For 
both  wavelengths  shown,  figure  6  can  be  interpreted  as  suggesting  that 
clearing  (associated  with  droplet  evaporation)  of  a  smoke  cloud  under 
conditions  of  decreasing  relative  humidity  is  slower  than  its  formation  at  any 
given  relative  humidity  between  about  50  and  75  percent.  The  figure  also 
demonstrates  that  the  sense  of  changes  in  relative  humidity  may  be  as 
important  as  its  value  in  modeling  hygroscopic  aerosols. 

In  figure  7  one  will  find  graphs  of  the  ratios  of  the  extinction  coefficients 
at  1.06pm  and  3.56pm  to  those  at  0.55pm  as  a  function  of  relative  humidity, 
figure  7  is  subject  to  various  interpretations.  One  inference  is  that  the 
aerosol  model  used  for  these  calculations  is  almost  always  (all  relative 
humidities  =  99  percent)  more  "transparent"  at  3.56pm  than  at  0.55pm,  but  the 
same  statement  cannot  be  made  at  1.06pm.  This  particular  aerosol,  while  more 
transparent  at  10.6pm  than  at  0.55pm  for  relative  humidities  below  about  59 
percent,  is  comparatively  less  transparent  for  larger  saturation  ratios. 
Other  features  discernable  from  figure  7  are  that  extinction  at  1.06pm  as 
compared  to  its  value  at  0.55pm  is  much  less  sensitive  to  relative  humidity 
than  it  is  for  3.56pm  radiation.  In  other  words,  the  extinction  at  1.06pm  can 
be  predicted  more  accurately  from  knowledge  of  the  extinction  at  0.55pm  than 
can  be  the  extinction  at  3.56pm  when  no  information  on  relative  humidity  is 
available.  One  further  obvious  inference  from  figure  7  is  that  extinction  at 
1.06pm  will  be  greater  than  at  0.55pm  if  the  smoke  cloud  is  clearing  because 
of  decreasing  relative  humidity. 


62 


Figure  7.  Ratio  of  Model  A'  extinction  coefficients  at  X  =  1.06ym 
and  X  =  3.56pm  to  values  at  X  =  0. 55Tiin  as  a  function  of 
relative  humidity.  The  solid  curves  represent  increasing 
relative  humidity,  and  the  dashed  curves  represent  de¬ 
creasing  relative  humidity. 
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Figure  8.  Attenuation  coefficient  per  unit  mass  of  Wet  Aerosol  at 
0.55ym  versus  Relative  Humidity  for  smoke  model  A'.  The 
solid  curve  is  for  increasing  relative  humidity  and  the 
dashed  curves  represent  decreasing  relative  humidity. 


Finally,  figure  8  shows  the  attenuation  coefficient  per  unit  mass  of  WET 
aerosol  at  X  =  0.55)im  as  a  function  of  relative  humidity.  If  that  figure  is 
compared  with  figure  6,  it  will  be  seen  that  the  attenuation  per  unit  mass  of 
the  existing  aerosol  (WET  -*■  nuclei  +  accreted  water)  is  much  less  sensitive  to 
relative  humidity  than  the  attenuation  per  unit  mass  of  dry  aerosol 
material.  It  suggests  that  simultaneous  in  situ  measurements  of  aerosol  mass 
concentrations  and  extinction  which  preserve  accreted  water  should  yield  less 
uncertainty  or  variation  in  extinction  with  changes  in  relative  humidity  than 
measurements  made  on  dry  aerosol  materials  alone. 

Discussion :  The  results  described  above  can,  of  course,  be  applied  in  detail 
only  to  the  aerosol  model  used  in  generating  them  and  at  the  specific 
wavelengths  used.  They  do  suggest,  however,  that  extinction  can,  under  some 
situations  at  least,  vary  by  nearly  a  factor  of  two  at  a  given  relative 
humidity  between  conditions  or  rising  and  falling  relative  humidity,  and  that 
the  relationship  between  extinction  at  two  different  wavelengths  may  also  be 
dependent  on  the  direction  of  changes  in  relative  humidity.  Consequently,  it 
would  seem  that  persons  involved  in  both  field  measurements  and  numerical 
modeling  should  pay  considerable  attention  to  meteorological  conditions,  in 
particular  when  the  relative  humidity  is  rising  or  falling. 

FURTHER  DESCRIPTION  OF  THE  "HALVING"  METHODS  OF  AGAUSX 

Basic  Concept.  A  general  idea  of  how  program  AGAUSX  operates  can  be  obtained 
by  considering  a  numerical  method  for  determining  the  area  under  a  curve  g(R) 
versus  R  such  as  that  sketched  in  figure  9.  The  objective  is  to  evaluate 
numerically  an  integral 


R=°° 

G  =  /  g(R)  dR 
R=o 


(61) 


to  some  desired  degree  of  accuracy  using  the  smallest  number  of  values  of  R 
for  the  numerical  calculations.  The  procedure  adopted  in  AGAUSX  is  as  follows 
(reference  to  figure  9  may  be  helpful): 
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(1)  An  initial  estimate  of  the  value  of  the  integral  is  made 
using  three  values  of  R  labeled  by  Roman  numeral  I  and  the  "trapezoidal  rule". 

(2)  A  second  estimate  G2  is  then  made  using  increments  AR  which  are 
half  as  large  as  those  used  in  getting  G^.  In  getting  G2>  the  two  additional 
R-values  labeled  II  are  utilized. 

(3)  The  values  of  G2  and  Gj^  are  compared  to  each  other  by 
calculating  a  quantity 


6  = 


(62) 


and  comparing  it  to  a  pre-set  quantity  A. 

If  6  <  A,  then  it  is  assumed  that  G2  is  a  "sufficiently"  accurate 
representation  of  G,  and  the  computations  are  terminated.  If  on  the  other 
hand,  6  >  A,  one  proceeds. 

(4)  A  third  estimate  Gg  of  G  is  then  made  by  again  cutting  the 
increments  AR  to  half  its  previous  value.  This  results  in  the  addition  of 
computations  at  the  four  new  R-values  labeled  III  in  figure  9.  A  new  value 
of  6  is  then  calculated  from 


6  = 


iG3-G2l 


(63) 


67 


and  6  is  again  compared  to  A .  If  6  is  still  greater  than  A,  the  spacing 
between  R- values  is  cut  in  half  once  more,  and  the  "estimation- comparison" 
process  is  repeated  until  either  (a)  6  <  A,  or  (b)  some  maximum  number  of  R- 
values  has  been  reached. 

In  AGAUSX,  the  quantity  used  in  the  testing  process  is  the  total  volume  of  the 
aerosol  particles,  namely, 


v  =  /  J  7TR3  f (R)dR, 


(64) 


in  which  f(R)dR  is  the  relative  number  of  aerosol  particles  whose  radii  lie 
between  R  and  R+dR. 

Extensions  to  the  Basic  Concept:  It  may  happen  that  the  function  g(R) 
illustrated  in  figure  6  will  have  a  form  which  is  much  "smoother"  over  some 
range  of  R-values  than  over  other  ranges.  Strict  application  of  the  halving 
procedure  to  such  a  situation  could  easily  lead  to  the  use  of  many  more  R— 
values  in  one  of  those  ranges  than  are  really  needed  because  of  behavior  of 
g(R)  in  other  R-ranges.  Such  a  situation  seems  especially  likely  with 
asymmetric  size— distributions.  One  way  of  attempting  to  avoid  unneeded 
computations  in  various  R-ranges  is  to  split  the  regions  used  in  halving  loops 
into  several  distinctive  "intervals"— RA  to  Rg,  Rg  to  Rc,  Rc  to  Rg,  and  to  do 
the  halving  operations  and  convergence  tests  separately  for  each  "interval." 
AGAUSX  has  been  coded  to  permit  separate  halving  computations  and  convergence 
testing  for  up  to  eight  "intervals"  of  that  type.  In  fact,  the  present  coding 
uses  two  radius  "intervals"  (ranges)  for  those  distributions  f(R)  known  to 
show  a  single  maximum,  four  "intervals"  for  bi-modal  distribution  models,  and 
just  one  interval  for  other  models. 
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DETAILED  DESCRIPTION  OF  PROGRAM  AGAUSX 


Program  AGAUSX  is  designed  to  calculate  various  light-scattering  quantities 
such  as  phase  functions  (Mie  theory),  scattering  fractions  (ACT),  extinction 
cross-section,  attenuation  coefficient,  etc. ,  for  diverse  natural  and 
artifically  created  polydisperse  atmospheric  aerosols.  The  program  consists 
of  subroutines  ANGLE,  GUSET,  AGXPT1,  AGXPT2,  AGXPT3,  AGXPRT,  MIEGS,  GAUS, 
VERFY,  WATER,  GPHASX,  GUSSTX,  and  PRINTX.  The  organization  and  operation  of 
these  subroutines  is  controlled  by  the  MAIN  program  AGAUSX.  Additionally 
program  AGAUSX  will  produce  an  analytic  (GHG)  phase  function  if  the  user 
desires. 


MAIN  PROGRAM  AGAUSX 

The  first  input  card  contains  the  parameters  NWAVE,  NINDX,  IDSTP,  NRADI,  IT, 
IANG,  NCRDS,  ITOT,  NUNIT,  MQRTE,  and  IAPX.  These  parameters  will  be  explained 
in  detail  subsequently.  The  various  modes  of  operation  of  AGAUSX  that  are 
possible,  are  controlled  by  some  of  these  parameters.  If  IANG  =  0,  the 
subroutine  GUSET  is  called  to  choose  'IT'  scattering  angles  between  0°  and 
180°  for  use  in  numerical  integrations  using  Gauss-Legendre  quadrature.  GUSET 
returns  two  arrays  of  numbers  to  MAIN.  The  array  C(I)  corresponds  to  the 
cosine  of  scattering  angles  and  the  array  H(I)  corresponds  to  the  quadrature 
we ight  s . 


If  IANG  =  1  or  2,  all  the  calculations  of  the  previous  paragraph  are 
skipped.  Instead  subroutine  ANGLE  is  called  and  returns  arrays  C(I)  and  H(I) 
to  MAIN.  For  IANG  =  1,  H(I)  corresponds  to  'IT'  equally  spaced  angles  between 
0°  and  180°,  and  C(I)  corresponds  to  the  cosine  of  these  angles.  If  IANG  =  2, 
'IT'  user  supplied  angles  are  read  in,  again  with  H(I)  corresponding  to  'IT' 
angles  and  C(I)  their  cosines. 

If  IAPX  is  greater  than  zero  AGAUSX  will  operate  in  a  special  'approximation' 
mode.  That  is  AGAUSX  will  construct  the  GHG  analytic  phase  function  at  IAPX 
Gauss-Legendre  angles.  When  approximate,  quick  phase  functions  are  desired, 
this  mode  should  be  used. 
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Next,  subroutine  AGXPT1  is  called.  One  or  more  input  parameters  are  read  at 
this  point  depending  on  the  value  of  IDSTP,  which  indicates  the  type  of 
distribution  to  be  used.  AGXPT1  calculates  and  returns  two  arrays  R(I)  and 
F(I)  which  describe  the  normalized  size  distribution,  and  VOL,  the  average 
'dry'  volume  per  particle.  The  array  R(I)  contains  values  of  'NRADI'  particle 
radii,  and  the  array  F(I)  contains  values  of  the  distribution  for  the 
corresponding  R(I).  For  IDSTP  =3,  7,  8,  9,  10,  11,  or  12  AGXPT1  also  returns 
DENS,  the  particle  number  density  (per  cm^).  For  other  distribution  types 
that  do  not  have  predetermined  density  values,  the  user  supplied  density, 
DENSH,  is  less  than  10""^;  in  this  case  the  particle  number  density  will  be 
calculated  from  mass  densities  and  mass  concentrations. 

The  parameters  WAVE,  EWAVE,  RELHUM,  DENSH,  TEMP  and  DELTA  are  read  next.  The 
next  do-loop  indexed  by  NWAVE  allows  computations  at  several  wavelengths 
during  a  single  run.  An  additional  looping  option  is  available  by  assigning 
DWAVE  any  value  less  than  10-^;  the  same  do-loop  is  used  to  handle  several 
values  or  relative  humidity  (RELHUM)  in  a  run. 

A  switch  parameter  LLLL  is  assigned  a  value  0  or  1,  in  a  manner  dependent  upon 
the  value  of  IDSTP  and  DENSH.  That  is  if  the  chosen  distribution  has  fixed 
parameters  and/or  if  number  densities  are  to  be  calculated  later  from  the 
average  particle  volume,  mass  density  and  mass  concentration,  then  LLLL  =  1. 

In  subroutine  AGXPT2  optical  and  physical  data  (index  of  refraction  and  mass 
density)  for  the  aerosol  material  are  read.  For  the  user's  convenience, 
subroutine  WATER  has  been  added  to  provide  the  relevant  data  for  water  limited 
to  the  range  0.2  to  200  micrometers.  If  RELHUM  and  the  growth  factor  (EMUA) 
are  non-zero,  AGXPT2  takes  into  account  the  particle  growth  and  the  change  in 
density  and  indices  of  refraction  of  the  aerosol  material  due  to  the 
absorption  of  water. 

In  AGXPT2,  subroutine  MIEGX  is  called  for  each  'adjusted'  radius.  MIEGX 
returns  to  AGXPT2  the  values  of  extinction,  scattering  and  backscattering 
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efficiency  factors,  and  the  average  intensity  factors  (1^  +  l2)/2  for  the  'IT' 
scattering  angles  chosen  earlier.  These  scattering  functions  are  weighted  and 
integrated  over  the  size  distribution. 

If  the  aerosol  is  a  mixture  of  components  having  different  physical  and/or 
optical  properties,  the  above  calculation  is  repeated  for  each  component  and 
various  functions  are  summed  over  NINDX  (the  number  of  components).  In  the 
end  AGXPT2  returns  to  MAIN,  for  each  wavelength,  the  average  (sum/total 
particle  number  density)  extinction,  scattering  and  back-scattering  cross- 
sections  (CTSUM,  CSSUM,  CRSUM),  the  average  intensities  (P(J)),  and  the  total 
mass  concentration  (TMASS,  in  gm/cc,  which  includes  the  mass  of  any  liquid 
water  absorbed  by  hygroscopic  aerosols). 

The  next  subroutine  is  called  AGXPT3.  AGXPT3  uses  the  numbers  received  from 
AGXPT2  and  returns  to  MAIN  the  total  extinction,  scattering,  and  back- 
scattering  coefficients,  scattering  fractions  and  phase  functions.  AGXPT3 
also  calculates  the  albedo  for  single  scattering  "ALBDO  “o  ^sca^ext  and 
prints  out  all  the  single  wavelength  results.  Later  in  AGXPT3,  If  IANG  =  0, 
subroutine  GAUS  is  called.  Using  the  quadrature  weights  H(I)  calculated 
earlier  and  calculating  the  Legendre  polynomials  PL(,),  GAUS  generates  the 
Legendre  expansion  coefficients  for  the  phase  functions,  and  places  them  in 
array  0L(I).  The  ai^'s  are  then  used  to  reconstruct  phase  functions,  PC(I): 
GAUSS  computes  the  root  mean  square  deviation  between  the  original  phase 
functions  and  the  reconstructed  phase  functions.  If  IANG  =  1,  the  above 
calculation  is  skipped. 

If  NWAVE  >  1,  all  the  calculations  of  AGXPT2  and  AGXPT3  are  repeated  (NWAVE-1) 
additional  times.  In  MAIN,  various  scattering  functions  are  summed  and 
divided  by  NWAVE,  and  subroutine  AGXPRT  is  called  to  print  out  all  the 
averaged  results.  AGXPRT  also  calls  GAUS  once  more  (if  IANG  =  0)  to  generate 
the  coefficients  for  the  averaged  phase  functions. 

If  IAPX  is  greater  than  zero,  subroutine  GPHASX  is  called  to  construct  an 
analytic  phase  function,  hereafter  referred  to  as  the  GHG  phase  function. 
This  requires  the  additional  subroutines,  GPHASX,  GUSSTX,  and  PRINTX.  This 
routine  will  construct  the  GHG  analytic  phase  function  at  IAPX  Gauss-Legendre 
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angles,  and  the  subsequent  Legendre  coefficients  may  be  written/punched 
unit  specified  by  NUNIT  if  the  user  so  desires. 


on  a 


SYMBOL 

ALBDO 

AM  AX 

C(  ) 

CATTN 

CRSUM 

CSSUM 

CTSUM 

DELTA 

DENS 


Explanation  of  Symbols  used  in  AGAUSX  (MAIN) 

_ Explanation  or  Definition _ 

single  scattering  albedo. 

the  largest  Mie  size  found  in  the  aerosol  distribution, 
array  holding  cosines. 

the  average  (sun/NWAVE)  attenuation  coefficient  in  square 
meters  per  milligram  of  aerosol  material. 

the  back-scattering  cross-section  in  square  micrometers  for 
each  wavelength  and  integrated  over  size  distribution. 
AGXPT3  returns  ! coefficients! . 

the  scattering  cross-section  in  square  micrometers  for  each 
wavelength  and  integrated  over  size  distribution.  AGXPT3 
returns  ’coefficients'. 

the  extinction  cross-section  in  square  micrometers  for  each 
wavelength  and  integrated  over  size  distribution.  AGXPT3 
returns  'coefficients'. 

convergence  criterion. 

the  particle  number  density  (number  per  cubic  centimeter) . 
The  value  of  DENS  may  be  supplied  by  the  user  as  an  input 
parameter,  DENSH.  However,  it  will  be  ignored  when  IDSTP 
=  0  or  greater  than  6,  because  these  distribution  types 
carry  predetermined  values  of  DENS.  In  the  case  of  IDSTP 
-  1,  2,  4,  5,  the  user  supplied  value  of  number  density 
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DENSH 

DRY  VOL 

DWAVE 

ELWC 

EMM 

ENWAV 

FSUM 

GNU 

H(  ) 

IANG 


will  be  ignored  only  if  it  less  than  10“^.  In  such  cases 
DENS  is  calculated  from  mass  density  and  concentration 
and  average  volume  per  particle,  and  is  represented  by 
DENSC. 

the  user-supplied  particle  number  density;  units  are 
particles  per  cubic  centimeter.  See  discussion  of  DENS 
(above) . 

the  average  volume  per  particle  of  dry  aerosol  in  cubic 
micrometers. 

the  wavelength  increment  in  micrometers,  or  relative 
humidity  looping  indicator. 

the  liquid  water  content  in  gm/cm^  (used  only  for  cases 
IDS TP  =  6  and  12). 

the  refractive  index  of  the  surrounding  medium. 

DFLOAT  (NWAVE) . 

the  numerical  result  for  the  integral  over  the  size 
distribution  function  with  respect  to  radius;  it  is  used 
to  normalize  the  distribution  function  to  one  (equivalent) 
particle  per  cubic  centimeter. 

wave  number  in  cm-^ . 

array  containing  angles. 

=  0,  for  the  computation  of  phase  functions  at  'IT1 
Gauss -Legendre  quadrature  angles. 
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IANG 

IDS  TP 

IAPX 

IDBLE 

IT 

ITOT 


KBAKT 

KEXTT 

KSCAT 


=  1,  for  computation  of  phase  functions  and  scattering 
fractions  at  l IT  *  equally  spaced  angles  between  0°  and 
180°  .  =  2,  for  computation  of  phase  functions  and 
scattering  fractions  at  'IT?  user  supplied  angles. 

identified  the  type  of  aerosol  size  distribution  to  be 
used.  It  can  only  take  values  between  0  and  12.  See 
AGXPT1  for  more  details. 

the  order  of  Legendre  expression  for;  the  analytic,  GHG, 
phase  function  expansion  when  IAPX  is  greater  than  zero. 

single  or  double  precision  mode  industor. 

the  order  of  Legendre  expansion  for  phase  functions  when 
IANG  =  0,  or  the  number  of  equally  spaced  angles  between 
0°  and  180°  when  IANG  =  1,  or  the  number  of  user  supplied 
angles  if  IANG  =  2. 

works  in  conjunction  with  NCRDS  and  NWAVE  to  write/punch 
either  individual  wavelength  values  of  the  phase  function, 
(Legendre  coefficients  and/or  scattering  fractions)  or 
averaged  wavelength  values  or  both  on  NUNIT:  ITOT  =  1 
for  individual  wavelengths;  ITOT  =  2  for  averaged  wave¬ 
lengths,  ITOT  =  3  for  both. 

the  average  (sum/NWAVE)  back-scattering  coefficient  per 
km,  integrated  over  the  size  distribution. 

the  average  (sum/NWAVE)  extinction  coefficient  per  km, 
integrated  over  the  size  distribution. 

the  average  (sum/NWAVE)  scattering  coefficient  per  km, 
integrated  over  the  size  distribution. 
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LLLL 


LMAX 


MQRTE 


NCRDS 


NINDX 


NRADI 


NUN  IT 


NWAVE 


OL(I) 


a  switching  parameter  used  to  control  whether  or  not 
particle  number  density  is  to  be  calculated  from  DRYVOL 
and  mass  concentration. 

=  3*IFIX  (AMAX).  Integer  estimate  of  the  optimal  order  for 
Gauss -Legendre  quadrature;  used  only  for  diagnostic  message 
print. 

=  12345,  will  cause  subroutine  AGXPT2  to  print  efficiency 
factors  QT,  QS,  and  QR,  and  the  normalized  distribution 
function  for  every  value  of  radius  used. 

=  1  for  write/punch  only  Legendre  coefficients  on  NUNIT : 

=  2  write  only  phase  functions  and  scattering  fractions  on 
NUNIT:  =  3  write/punch  both  on  NUNIT. 

the  number  of  aerosol  components  which  will  have  different 
optical  constants,  mass  density  or  mass  concentration. 

the  number  of  values  of  particle  radius  to  be  used  in  the 
calculations.  (In  effect,  points  on  the  radius  versus  size 
distribution  function  plot). 

defines  the  device  on  which  output  data  may  be  stored  or 
punched;  may  be  used  to  place  nominal  card  output  into 
data  files,  on  tape,  etc.  The  default  value  (NUNIT  =  0) 
is  4. 

the  number  of  wavelengths  or  relative  hunidity  values  to  be 
treated  in  a  given  run.  EWAVE  has  to  be  less  than  10-^  for 
the  latter. 

the  average  Legendre  expansion  coefficients. 

OLT(I)/ENWAVE. 
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OLT(I) 


the  total  Legendre  expansion  coefficients  (summed  over  all 
values  of  !WAVE!). 


OUT(  ) 


an  array  used  for  storing  some  of  the  numbers  to  be  printed 
later. 


PL(I,J) 


the  Legendre  polynomials  of  order  (1-1)  and  argument  C(J). 
C(I)  are  cosines  of  the  scattering  angles. 


PSUM(I) 


the  average  phase  functions  integrated  over  the  size 
distribution. 


PSUMT(I) 

QATTEN 

R(NRADI) 

RELHUM 

SCAT(I) 

SCATT(I) 

TEMP 

TMASS 

VOL 

WAVE 


the  total  phase  functions  integrated  over  size  distribution 
and  summed  over  wavelength. 

the  attenuation  coefficient  in  square  meters  per  milligram 
of  aerosol  material  for  each  wavelength. 

the  radius  of  the  largest  particle  encountered. 

the  relative  humidity  in  percent. 

the  scattering  fractions  as  defined  in  ACT  for  each 
wavelength.  'IT!  elements  in  the  array. 

the  average  (sum/NWAVE)  scattering  fractions  (ACT). 

the  temperature  of  the  atmosphere  in  degrees  C. 

the  total  mass  of  aerosol  in  gm/cm^  in  the  atmosphere. 

the  average  volume  per  particle  in  cubic  micrometers. 

the  first  wavelength,  in  micrometers,  at  which 
calculations  are  to  be  done. 
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AG AUS X  (MAIN  PROGRAM)  -  SIMPLIFIED  FLOWCHART 
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SUBROUTINE  AGXPT1 


Readers  of  the  following  descriptions  of  the  mathematical  procedures  should 
bear  in  mind  that  the  general  AGAUS  single  scattering  code  computes  a  number 
of  quantities  associated  with  the  Mie  theory  of  the  interaction  of  a  single 
spherical  scatter  and  then  averages  those  quantities  using  the  aerosol  size 
distribution  function,  f(r),  as  a  weighting  factor.  The  following  definitions 
may  also  be  helpful. 

The  terms  "component"  or  "aerosol  component"  apply  to  cases  in  which  the 
overall  aerosol  may  be  a  mixture  of  different  materials  whose  optical 
properties,  mass  density,  growth  factors,  or  mass  concentrations  are  not 
identical.  The  number  of  "components"  to  be  used  in  a  run  is  specified  by 
input  parameter  NINDX.  An  example  of  a  multi-component  aerosol  might  be  one 
which  is  a  mixture  of  hydroscopic  and  nonhydro scopic  particles. 

The  terms  "size  interval"  and/or  "original  interval”  refer  to  subdivisions  of 

the  total  range  of  aerosol  radii  (Vnimum  =  RL0  to  Maximum  =  RHI) 

which  the  halving  and  convergence  testing  procedures  are  applied 

independently. 

Most  averages  that  must  be  done  are  of  the  form 


RHI 

G  =  J  dr  f (r)G(r)/FSUM,  (67) 

RLO 


where  G(r)  is  some  Mie  theory  quantity  for  a  particle  of  radius  r,  and  a  given 
optical  type  (refractive  index  m  =  m-ik) ,  and 


RHI 

FSUM  =  /  dr  f (r) . 

RLO 
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The  size  distributions  f(r)  ,  RLO  <_  r  RHI,  may  be  given  in  analytic  form,  or 
numerically  on  a  set  {R(J)}  of  values  of  r.  They  may  or  may  not  be  normalized 
to  unity  (FSUM  may  equal  1.0,  or  not).  The  relevant  Mie  theory  quantities 
include  cext(r)>  Csca^r^»  P(r>v)  for  each  Gauss -Legendre  p,  etc.  Other 
averages  that  must  be  done  have  somewhat  different  forms,  e.g., 


ft) 

n 


dr  f(r)wn(r)Cext(r)/FSUM, 


(68) 


for  the  coefficients  co  of  the  phase  function  expression.  p(y)  = 

I  terms  of  the  coefficients  oo  (r)  of  the  phase  function 

n=o  n 

00 

p(r,y)  =  l  wn(r)Pn(y). 


The  above  integrals  over  r  must  be  done  numerically.  In  order  to  minimize 
computation  time,  the  number  of  points  in  the  numerical  integration  should  be 
chosen  as  small  as  possible,  consistent  with  adequate  accuracy.  Each  value  of 
r  that  is  used  requires  the  full  Mie  calculation  for  the  quantities  C  (.(r)  , 
Cer.Q(r)>  p(r>u)  for  each  p,  etc.,  these  calculations  absorb  a  large  part  of 
the  total  computation  time.  It  has  been  found  in  this  work  that,  instead  of 
fixing  the  number  of  values  of  r  to  be  500,  adequate  accuracy  results,  in  many 
cases,  for  as  few  as  50  to  100  values  of  r;  this  reduction  allows  a  reduction 
in  total  computation  time  by  roughly  a  factor  of  4  to  8. 

The  numerical  integration  method  which  was  developed  makes  use  of  successive 
halving  of  intervals  until  a  preset  convergence  criterion  is  met.  This 
halving  method  allows  an  initial  set  of  unequal  intervals  to  be  chosen.  For 
example,  suppose  that  the  distribution  function  f(r)  is  sharply  peaked  around 
some  value  of  r,  and  also  has  a  long  tail,  as  in  figure  10.  Hien  it  makes 
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RR(1) 

RLO 


RR(2) 


RR(3)  r 
RHI 


FIGURE  10;  Sketch  of  a  hypothetical  size-distribution. 

The  logical  flow  through  the  halving  integration  of  AGXPT1  is  as 
follows : 

(1)  Given  RLO,  RHI  and  f(r),  RLO  £  r  £  RHI,  in  analytic  form;  or 
f(r)  numerically  on  the  equally  spaced  set  r  -  R(J) ,  J  -  1,  NRADI. 

(2)  For  the  given  f(r),  RLO,  RHI,  choose  the  number  and  size  of  the 
original  intervals.  That  is,  choos  RR(I) ,  1=1,  NLAST  =  2**MIN  +  1, 
with  RR(1)  =  RLO,  RR(NLAST)=  RHI.  Note  that  NI  =  number  of  initial 
intervals  =  2**MIN  =  1,2,4,8,16,...  .  Then  set  NHALV  =  number  of  halvings 
in  each  interval  =  MAX-MIN,  and  NK  =  maximum  number  of  points  in  each 
interval  =  2**NHALV. 

(3)  The  maximum  number  of  points  allowed  in  NMAX  =  1+2**MAX.  In 
this  work,  MAX  =  9  was  used,  so  NMAX  =  513.  Calculate  R(KK) ,  F(KK)  = 
f(R(KK)),  KK  =  1,  NMAX,  on  all  the  points  R(KK)  which  might  be  used  if 


the  halving  goes  all  the  way  to  NMAX  points.  For  f(r)  given  numerically, 

this  is  done  by  linear  interpolation. 

RRfI+1) 

(4)  Calculate  C  (I)  =  J  dr  f(r)C  (r)  in  each  interval  number 

RR(I)  ext 

I,  I  =  1,  NI,  by  halving  of  that  interval  and  the  trapezoidal  rule,  until 
the  convergence  criterion  is  met  in  that  interval,  or  until  the  halving  has 
progressed  all  the  way  to  NK  points,  the  maximum  number  allowed  in  any  interval. 


RRfI+1) 

At  the  same  time,  calculate  FSUMG(I)  =  J  dr  f  (r)  ,  and  all  other  needed 

RR(I) 


RRfI+1) 

averages  G(I)  =  J  dr  f(r)G(r),  on  the  same  set  of  points  used  for  C  (I). 

RR  Cl)  6Xt 


(This  means  that  the  convergence  criterion  is  applied  only  to  C  (I)).  Then 

6  xt 

NI  NI 

Cext  =  l  Cext(I)>  FSUM  =  l  FSUM(I) ,  etc.;  then  Cext  =  C^/FSUMG, 

RHI 

G  =  G/FSUMG.  The  last  step  is  equivalent  to  making  J  dr  f(r)  =  1.0. 

RLO 
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sense  to  choose  unequal  initial  intervals  as  shown.  The  numerical  integration 
then  proceeds  by  halving  each  of  the  original  intervals,  until  the  convergence 
criterion  is  met.  The  convergence  criterion  that  was  used  is  the  requirement 
that  |G^n+^  -  G^n)|/|G^n^|  <  A,  where  A  is  some  small  preset  number.  Here 
is  the  value  of  the  integral  G  after  the  nth  halving,  separately  in  each 
initial  interval.  The  trapezoidal  rule  was  used  throughout.  The  value  of  the 
integrand  at  any  given  point  r  is  calculated  just  once;  it  does  not  have  to  be 
recalculated  on  that  point  after  halving. 

Description  of  Types  of  Distributions 

IDSTP  DESCRIPTION 

0  This  is  an  arbitrary  user-supplied  distribution.  'NRADI'  +  1 

cards  will  be  read;  the  first  card  contains  RLO  and  DELLR  (vim). 

The  rest  of  the  cards  carry  the  values  of  F(J)  and  must  be  in 
order  of  increasing  radius  value. 

1  This  is  the  zero-order  log-normal  distribution:  the 

distribution  function  is  given  by 


F(R)  = 


Hog  (R/R)  „ 


/2tT  Hoge(a)R 


eXI><  -  2  hog »  1  1 


R  =  RBAR;  a  =  SIGMA,  is  the  standard  deviation.  This 
distribution  type  requires  one  input  data  card  to  read  in  the 
values  of  R,  0,  RLO  and  RHI. 
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2  This  is  called  the  double  exponential  distribution  and  its 
distribution  function  is  given  by  F(R)  =  QA  exp(-AR)  + 

(l-Q)B  exp(-BR).  Q  5  CUE.  This  distribution  type  requires 
one  input  data  card  to  read  in  the  values  of  RLO,  RHI,  Q,  A, 
and  B.  Q  is  dimensionless  while  A  and  B  have  units  of 

3  This  model  (Deirmend jian' s  "Model  C")  does  not  require  any 
input  data  card.  It  carries  fixed  value  of  DENS,  RLO,  and 
DELRD.  RHI  is  determined  by  the  input  parameter  NRADI. 

F(R)  =  450.2  R  <  .08 

=  2.251*DELRD*R~4  R  >  .08 

4  The  distribution  function  of  this  model  (Junge  distribution) 
is  given  by  F(R)  =  QR~\  Q  =  CUE.  This  distribution  type 
requires  one  input  data  card  to  read  in  RLO,  RHI,  A,  A. 

5  The  distribution  function  for  the  Modified  Gamma/ Generali zed 
Khirgian-Mazin  distribution  is 

F (R)  =  Ra  exp  ,  a  =  ALF,  y  =  GAM,  and 

c  ' 

Rc  s  RC.  One  input  card  is  needed  to  read  in  RLO,  RHI,  RC, 

ALF,  GAM,  and  ELWC.  EIWC  is  not  needed  for  type  5  distribution 
and  therefore  can  be  left  blank. 

6  The  size  distribution  model  (NMSU  Fog  or  Cloud  Model)  is  very 
similar  to  type  5,  except  that  the  user  must  supply  one 
additional  input  parameter — namely,  the  liquid  water  content 
(EIWC)  in  gm/cc.  This  model  can  be  used  for  treating 
situations  involving  liquid  water  aerosols  like  clouds  or  fogs. 
For  type  6  runs  one  does  not  need  to  read  in  the  values  of  EMA, 
CAYA,  RHOA,  CONC. 
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This  distribution  is  essentially  same  as  Junge's  distribution 
(type  4)  except  that  it  has  fixed  parameters.  One  input  card 
is  needed  to  read  in  VIS  (visibility  in  km);  VIS  is  used  in 
calculating  DENS. 

This  is  a  fixed  parameter  Continental  Bi -modal  Model.  It  does 
not  require  an  input  data  card. 

This  is  a  fixed  parameter  Maritime  Bi -modal  Model.  It  does  not 
require  an  input  data  card. 

This  is  a  fixed  parameter  Urban  Bi -modal  Model.  It  does  not 
require  an  input  data  card. 

This  is  a  user-supplied  Bi -modal  Model.  This  requires  one 
input  card  to  read  in  FOA,  FOC,  SGA,  SGC,  RBARA,  RBARC.  Types 
8,  9,  10,  and  11  use  the  sum  of  two  log-normal  distributions: 

1  ^(R/R.)  2 

2  ^oge(ai) 

N^  =  FOA,  R^  =  RBARA,  =  SGA  with  similar  meaning  for  N2, 

R2,  and  02*  Note  that  in  Type  8,  9,  10  the  values  of  SGA 
and  SGC  are  £oge(a). 

This  model  (Marshall-Palmer  Rain  Model)  is  a  simple  exponential 
model  which  assumes  an  empirical  relationship  between  rain  rate 
and  droplet  size  distribution  parameters: 


F(R)  =  l 


N. 

1 


i=l  /2tt  iloge(ai)R 


exp 


SYMBOL 

A 

ALF 

AVOL 

DENA 

DENC 

DENS 

DELR 

DELRD 

DELLR 


NQ  =  0.08  cm  4,  and  A  =  41(RN)-0,21  in  which  RN  =  RAIN  is  the 
rain  rate  in  mm/hour.  Diameter  D  is  in  cm.  The  corresponding 
size  distribution  function  of  radius  R  is  given  by 


F (R)  =  2N  exp(-2AR) . 
o 


This  distribution  requires  one  input  data  card  to  read  in  RAIN. 
The  values  of  RLO  and  RHI  are  fixed  at  0  and  0. 5  cm 
respectively.  Due  to  the  limitations  on  the  range  of  Mie- 
sizes  (subroutine  MIEGX)  type  12  usage  is  limited  to  wave¬ 
lengths  of  the  order  of  1  mm  or  larger.  Since  subroutine 
WATER  does  not  contain  optical  data  for  wavelengths  longer 
than  0.2  mm,  type  12  runs  require  the  user  to  supply  the 
values  of  EMA,  CAYA,  and  RHOA  as  if  rain  were  a  non-aqueous 
aerosol. 


_ Explanation  or  Definition _ _ 

parameter  in  several  distribution  types. 

parameter  in  modified  gamma  distribution. 

the  average  volume  per  particle  in  cubic  microns  obtained 

via  analytical  integration  over  the  limits  RLO  =  0  and  RHI  =  00 . 

That  has  only  been  done  for  IDSTP  =  5,6,8,9,10,11,12. 

temporary  storage. 

temporary  storage. 

_3 

the  particle  number  density  in  cm 

increment  of  radius  for  arbitrary  distribution. 

=  (RHI-RLO) /RADS;  increment  between  successive  values  of  R. 
increment  in  radius  for  the  case  IDSTP  =  0. 
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SYMBOL  _ Explanation  or  Definition _ 

DR  (  )  initial  size  ingerval  I  =  1,  NI. 

EM  real  part  of  refractive  index  of  actual  aerosol  particle. 

EMM  refractive  index  of  medium  in  which  scattering  particles  are 

deployed. 

F(J)  the  array  containing  'NRADI '  values  of  the  size  distribution 

function.  See  the  description  of  AGXPT1  for  more  details. 

FKK  value  of  distribution  function  at  radius  R(KK)/ 

FSUM  the  numerical  integral  over  the  size  distribution  function 

with  respect  to  radius  between  the  limits  RHI  and  RLO;  used  to 
normalize  the  distribution  function. 

GAM  parameter  in  modified  gamma  distribution. 

GNUM  temporary  storage. 

GNUMA  temporary  storage. 

GNUMC  temporary  storage. 

IDSTP  identifies  the  type  of  aerosol  size  distribution  to  be  used. 

IW  switch  for  "all  water"  mode. 

NCRDS  not  used  in  this  subroutine. 

NHALV  maximum  number  of  interval  halvings  allowed  for  each  initial 
interval . 

NI  number  of  initial  intervals  for  halving  integration. 

NKG  maximum  number  of  points  which  may  be  used  in  each  initial 

interval . 

NLAST  number  of  radii,  RR Cl),  defining  the  NI  basic  size  intervals. 

NMAX  maximum  number  of  points  which  may  be  used  in  halving  method. 


*The  pre-coded  (IDSTP  =  8.9,10)  values  of  SGA  and  SGC  are  the  natural 
logrithms  of  the  standard  deviations. 
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SYMBOL 


NRADI 

RAIN 

RBAR 

RBARA 

RBARC 

RC 

RHI 

R(J) 

RLO 

RR(I) 

SGA 

SGC 

SIGMA 

VIS 

VOL 

VOLA 

VOLC 


_ Explanation  or  Definition _ 

the  number  of  radius  values  to  be  used  in  describing  the  size 
distribution  function  for  type  zero  and  3  only, 
parameter  in  Marshall -Palmer  rain  model, 
mean  radius  in  log-normal  distribution. 

mean  radius  of  "accumulation"  mode  for  bimodal  distribution, 
mean  radius  of  "coarse"  mode  for  bimodal  distributions, 
mode  radius  in  modified  gamma  distribution, 
maximum  particle  radius  for  any  size  distribution  (pm) . 
the  array  containing  'NRADI'  values  of  radius  (particle  size) 
in  micrometers. 

minimum  particle  radius  for  any  size  distribution  (pm) . 
lower  radius  of  Ith  initial  interval:  RR(1)  =  RLO,  RR(NLAST)  =  RHI. 
standard  deviation*  for  user  supplied  bi-modal  log  normal  distribution, 
standard  deviation*  for  user  supplied  bi-modal  log  normal  distribution. 


standard  deviation  in  log-normal  distribution:  used  later  as 
In  (SIGMA), 
visibility  in  km. 

3 

the  average  'dry'  volume  per  particle  (in  pm  )  calculated 

numerically. 

particle  volume. 

particle  volume. 
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Subroutine  AGXPT1  -  Simplified  Flow  Chart 


RETURN 
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SUBROUTINE  AGXPT2 


Subroutine  AGXPT2  deals  with  the  most  important  computational  aspects  of 
AGAUSX.  If  RELHUM  and  the  growth  factor  (EMUA)  are  non-zero,  AGXPT2  makes 
adjustments  in  particle  radii,  density  and  refractive  index  of  aerosol 
material  to  account  for  the  absorption  of  water  by  aerosol  particles.  These 
adjustments  are  made  through  the  following  equations: 

Radius:  RT  =  AC*R(L)  -  BC/AC, 

Index  of  Refraction: 

i)  Real  part:  EM  =  EMW  +  (EMA-EMW)/A, 
ii)  Imaginary  part:  part  CAY  =  CAYW  +  (CAYA  -  CAYW)/A 
Density:  RHO  =  RHOW  +  (RHOA  -  RHOW)/A, 


All  the  symbols  used  above  are  explained  in  detail  below. 

Next,  subroutine  MIEGX  is  called  to  compute  single  particle  scattering 
functions  QR,  QS,  QT,  and  P(J)  for  each  adjusted  radius.  These  scattering 
functions  are  weighted  and  integrated  over  the  size  distribution  by  the 
trapezoidal  method.  The  convergence  of  integration  for  volume  and  the 
extinction  cross-section  is  checked,  and  warnings  are  printed  if  the  final  one 
or  two  contributions  exceed  5  percent  of  the  previous  total  values.  If  there 
is  more  than  one  component  (NINDX>1)  in  the  aerosol  with  different  refractive 
index  and  density,  the  above  calculation  is  carried  out  NINDX  times,  and 
results  are  summed  and  divided  by  the  total  particle  number  density.  In  the 
end,  the  attenuation  coefficients  in  square  meters  per  milligram  of  'wet' 
aerosol  material  are  calculated.  The  same  value  of  extinction  coefficient  is 
used  for  calculating  both  attenuation  coefficients. 

If  IAPX  is  greater  than  zero,  only  the  first  three  coefficients  of  the 

Legendre  expansion  (wQ,  u^)  are  calculated  and  stored  in  array 

0L(  ).  This  is  for  subsequent  construction  of  the  GHG  analytic  phase 

function. 
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SYMBOL 


Definition 


A  =  1  +  (RHOA/RHOW)  *  EMUA  *  CH;  used  in  making  adjustments  due 

to  water  absorption  by  aerosol. 


AC 

=  A1/3. 

ALBDO 

single  scattering  albedo. 

ALPHA 

=  2*PI*EMM*RT/WAVE.  Mie  size  parameter 

for  adjusted  radii. 

ALPHAD 

temporary  storage  for  ALPHA. 

BC 

=  BHT*CH;  used  in  making  adjustments 

in 

size  growth. 

BH 

=  1.056*10“3. 

BHT 

=  BHI (298/TEMK) . 

CATTN 

the  total  attenuation  coefficient  in 

square  meters  per  milligram 

of  'dry'  aerosol  material  for  each  WAVE.  However,  CTSUM 
corresponds  to  'wet'  aerosol. 


CATTNW  the  total  attenaution  coefficient  in  square  meters  per  milligram 
of  'wet'  aerosol  material. 

CAY  the  imaginary  part  of  refractive  index  of  dry  aerosol. 

(Assumed  to  be  negative;  do  not  enter  a  value  with  a  negative 
sign  on  the  data  card(s)). 

CAYO  temporary  storage. 

CAYW  the  imaginary  part  of  refractive  index  of  pure  water  at  a  given 

TEMP  (DEG  C)  and  WAVE. 

CH  =  FH/l-FH) . 

CONC  the  mass  concentration  in  gm/cc  of  a  component  of  dry  aerosol. 

It  is  the  number  of  grams  of  dry  aerosol  per  cubic  centimeter 
of  "cloud"  or  "fog",  etc. 

CONCT  the  total  mass  concentration  in  mg/cc  of  dry  aerosol. 

CRGG  partial  contributions  to  average  radar  cross  section,  for  one 
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SYMBOL 


CRHH 

CRSUM 


CRSUMT 

CJGG, 

CSHH 

CSSUM 


CTSUM 

CTSUMT 

DEL 

DELTA 

DENS 

DENSC 


DENST 

DRYVOL 

ELWC 


_ Definition  _ 

component  of  aerosol. 

the  average  (sum/DENST)  back-scattering  cross-section  in 
square  micrometers,  integrated  over  the  size  distribution,  for 
each  WAVE. 

value  of  radar  corss  section  for  total  aerosol. 

partial  contributions  to  average  scattering  cross  section  for 

one  component  of  aerosol. 

the  averge  (sum/DENST)  scattering  cross-section  in  square 
micrometers,  integrated  over  the  size  distribution,  for 
each  WAVE. 

the  average  (sum/DENST)  extinction  cross-section  in  square 
micrometers,  integrated  over  the  size  distribution,  for  each  WAVE, 
the  total  extinction  cross-section  in  square  micrometers, 
integrated  over  the  size  distribution,  for  each  WAVE, 
fractional  change  in  contribution  to  volume  from  n^  to  n+lSt 
halving,  for  one  initial  interval, 
convergence  criterion. 

the  particle  number  density  per  cubic  centimeter.  See  DENSC 
also. 

the  particle  number  density  per  cc  for  each  component  of  aerosol. 
It  is  calculated  from  mass  density  and  concentration,  and  average 
volume  per  particle.  If  LLLL=1 ,  the  calculated  value  of  DENSC 
is  replaced  by  DENS,  the  value  of  which  has  been  determined  or 
supplied  elsewhere. 

the  total  particle  number  density  per  cc.  See  also  DENSC. 
the  average  volume  per  particle  in  cc. 
the  liquid  water  content  in  gm/cc. 
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SYMBOL 


EM 

EMA 

EMASS 

EMM 

EMUA 

EMW 

F  ( I ) 
FF  (J) 

eft 

FH 

FKK 

FKKA 

FT 

H(  ) 

IAPX 

IDBLE 

IDSTP 


_ Definition _ 

the  real  part  of  the  effective  refractive  index  of  'wet' 
aerosol.  See  CAY  EM  =  (EMW  +  (EMA  -  EMW) /A) . 
the  real  part  of  refractive  index  of  the  dry  aerosol, 
is  'wet'  mass  concentration  in  gm/cc  of  a  component  of  an 
aerosol  when  RELHUM  and  growth  factor  are  non-zero.  EMASS  and 
CONC  should  have  some  value  if  growth  factor  is  zero, 
the  refractive  index  of  the  surrounding  medium,  set  equal  to  1 
here  implying  no  surrounding  medium. 

Hanel's  mass  accretion  coefficient  y.  Values  must  ve  supplied 
by  the  user  and  depend  on  the  type  or  composition  of  the  aerosol 
being  modeled  as  well  as  upon  the  value  of  the  relative  humidity, 
the  real  part  of  refractive  index  of  water  at  a  given  TEMR 
(DEG  C)  and  WAVE. 

the  normalized  size  distribution  function  (in  micrometers-*), 
value  of  arbitrary  distribution  functions  at  radius  RP  = 

RLO  *  (J-l)  +  DELLR,  J  =1,  NRADI . 

contribution  to  integral  of  destribution  function  from  one 
original  interval. 

fractional  relative  humidity  (saturation  rates) . 
value  of  distribution  function  at  radius  R(KK) . 
geometrical  cross  section  of  particle. 

partial  contribution  to  integral  of  distribution  function, 
array  containing  angles. 

order  of  expansion  for  analytic  phase  function, 
precision  mode  indicator. 

identifies  the  type  of  aerosol  size  distribution  being  used  in 
the  run. 
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SYMBOL 


I ERROR 


KEXOLD 


KEXT 

KEXTT 

LLLL 


MERROR 

MQRTE 

NCRDS 

NHALV 

NINDX 

NRADI 

01 STAR 
01STRD 
02STAR 
02STRD 


_ Definition _ 

a  flag  which  is  set  to  unity  if  the  number  of  terms  reaches  the 

maximum  value  allowed  in  the  dimensions  of  the  Mie- 

coefficients  a  and  b  . 

n  n 

the  extinction  coefficient  per  km  summed  over  the  number  of 
components  one  less  than  the  current  value  of  the  running  loop 
index  NK. 

the  extinction  coefficient  per  km  for  each  component  of  the 
aerosol . 

the  total  extinction  coefficient  per  km  for  each  WAVE, 
a  switch  parameter.  If  LLLL=Q,  the  particle  number  density 
(DENSC)  is  calculated  using  mass  density  and  concentration,  and 
the  average  volume  per  particle  (DRYVOL) .  If  LLLL=1,  a  pre¬ 
calculated  or  pre-supplied  value  of  DENS  is  used, 
an  error  counter:  If  MERROR  exceeds  10,  execution  is  terminated. 

=  12345,  QT,  QS,  QR  and  Fr  are  printed  for  each  radius, 
not  used  in  this  subroutine,  but  appears  in  a  common  block, 
maximum  number  of  interval  halvings  allowed  for  each  initial 
interval . 

the  number  of  aerosol  components  which  will  have  different 
optical  constants,  mass  density  or  mass  concentration, 
the  number  of  points  on  the  radius  vs.  size  distribution 
function  plot. 

value  of  for  a  given  size  parameter, 
temporary  storage  for  01STAR. 
value  of  for  a  given  size  parameter, 
temporary  storage  for  02STAR. 
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SYMBOL 


Definition 


0L(  ) 

0L1GG, 

0L1HH 

OL1HHT 

OL1SUM 

OL2GG, 

OL2HH 

OL2HHT 

OL2SUM 

OLSTAR 

OM2 

P(J) 

PGG(I), 

PHH(I) 

PHHT(I) 

PL(  ,  ) 

PSNEW(J) 

PSOLD(J) 


PSUM(J) 

PSUMT(J) 

QR 

QRD 


Legendre  coefficients. 

partial  contributions  to  average  to  ,  for  one  component  of 
aerosol . 

contributes  to  to^  of  average  particle  from  one  original  interval, 
value  of  to  for  one  component  of  an  aerosol. 

partial  contributions  to  average  for  one  component  of  aerosol 

contributions  to  of  average  particle  from  one  original  interval, 
value  of  co £  for  one  component  of  aerosol . 
second  Legendre  coefficient  (co^) . 
third  Legendre  coefficient  (o^) . 

an  array  containing  'IT'  average  intensity  factors  +  i2)/2  for 
each  radius. 

partial  contributions  to  average  intensity  (i^  +  i^) /2  for  one 

aerosol  component  at  abscessa  value  p  . 

average  intensity  at  p^.  for  one  original  interval. 

Legendre  polynomials. 

th 

=  P(J)*F(L),  where  p(J)  corresponds  to  L  radius,  L  2. 

t  h 

=  P(J)*F(L),  where  P(J)  corresponds  to  (L-l)  radius  except 
when  L=1 ;  in  that  case  PSOLD(J)  =  P(J)*FF(1),  and  P(J) 
corresponds  to  1st  radius. 

average  intensity  factors  integrated  over  the  size  distribution. 
They  are  then  summed  over  NINDX  components  and  divided  by  DENST. 
final  average  phase  function. 

the  back- scattering  efficiency  factor  for  each  radius, 
temporary  storage  for  QR. 
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SYMBOL 


QS 

QSD 

QT 

QTD 

R(L) 

RELHUM 

RHO 

RHOA 

RHOW 

RIT 

RR(l) 

TEMK 

TEMP 

TMASS 

TVOL 

VOL 

VOLGG 

VOLHH 

VOLHHT 

WAVE 


_ Definition _ 

the  scattering  efficiency  factor  for  each  radius, 
temporary  storage  for  QS. 

extinction  efficiency  factor  for  a  given  size  parameter, 
temporary  storage  for  QT. 

the  array  containing  'NRADI'  values  of  radius  in  micrometers. 

the  relative  humidity  in  percent. 

the  specific  density  (in  gm/cc)  of  'wet'  aerosol. 

the  specific  density  (in  gm/cc)  of  'wet'  aerosol. 

the  specific  density  (in  gm/cc)  of  water  at  the  temperature 

TEMP  (in  DEG  C) . 

the  radius  of  a  particle  after  taking  into  account  its 
growth  due  to  absorption  of  water. 

lower  radius  of  I*"*1  initial  interval. 

the  temperature  of  surrounding  medium  in  degrees  Kelvin, 
the  temperature  of  the  surrounding  medium  in  degrees  Centigrade, 
the  total  mass  concentration  in  gm/cc  of  'wet'  aerosol, 
is  used  to  pass  the  value  of  DRYVOL  from  MAIN  to  subroutine. 

3 

total  volume  occupied  by  aerosol  material  distributed  in  1  cm 
of  space. 

partial  contributions  to  volume  of  average  particle,  for  one 
component  of  aerosol. 

contribution  to  volume  of  average  particle  from  one  original 
interval . 

the  wavelength  in  micrometers  at  which  all  the  scattering 
functions  are  computed. 


Subroutine  AGXPT2  -  Simplified  Flowchart 
(Receives  control  data  and  size  distribution  data  from  main  program) 
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SUBROUTINE  AGXPT3 


Subroutine  AGXPT3  receives  the  values  of  various  cross-sctions  and  the 
integrated  average  intensity  factors  from  AGXPT2  and  converts  them  into 
directly  useful  quantities: 


Coefficient  (per  km)  -  10  ^*DENS*cross-section, 

Scattering  fractions  =  — j  *DENS*intensity  factors, 

4tt 

X2 

and  Phase  function  =  — -  *intensity  factors. 

ext 


X  is  in  pm.  It  then  prints  out  all  single  wavelength  results. 


It  also  computes  albedo  for  single  scattering  using  equation  (4) .  It  then 
calls  on  subroutine  GAUS  which,  among  other  things  (see  GAUS  for  details), 
also  returns  a  value  for  albedo  calculated  differently.  The  two  values  are 
compared,  and  if  they  differ  by  more  than  0.01  percent  the  user  is  advised  to 
rerun  the  program  using  a  larger  value  of  IT. 


SYMBOL 

ALBEDO 

C(I) 

CAY 

CAYNG 

CRSUM 


CSSUM 


_ Explanation  or  Definition _ 

=  CSSUM/CTSUM;  albedo  for  single  scattering, 
the  cosines  of  scattering  angles.  'IT'  elements  in  the  array, 
the  ratio  of  the  imaginary  part  to  the  real  part  of  adjusted 
refractive  index  of  the  last  component. 

=  -CAY. 

the  average  back-scattering  cross-section  when  it  is  received  by 
AGXPT3  but  it  returns  to  MAIN  the  value  of  the  average  scattering 
coefficient  per  km. 

the  average  scattering  cross-section  when  it  is  received  by 
AGXPT3  but  it  returns  to  MAIN  the  value  of  the  average 
scattering  coefficient  per  km. 
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SYMBOL 

CTSUM 


DENS 

EM 

EMM 

GNU 

H(I) 

IANG 


IDSTP 

IT 


ITOT 

KRSUM 

KSSUM 


_ Explanation  or  Definition _ _ 

the  average  extinction  cross-section  when  it  is  received  by 
AGXPT3  but  it  returns  to  MAIN  the  value  of  the  average  extinction 
coefficient  per  km. 

_3 

the  total  particle  number  density  in  cm 

the  real  part  of  the  adjusted  refractive  index  of  the  last 
aerosol  component. 

the  refractive  index  of  the  surrounding  medium.  It  is  set  equal 

to  1  in  AGAUSX. 

the  wave  number  in  cm  * . 

the  array  containing  ’IT'  scattering  angles  (in  degrees). 

=  0,  for  computation  of  phase  functions  at  'IT'  Gauss- Legendre 
quadrature  angles. 

=  l,  for  computation  of  phase  functions  and  scattering 
functions  at  'IT'  equally  spaced  angles  between  0°  and  180°. 

=  2,  for  computation  of  phase  functions  and  scattering  fractions 
at  'IT'  user  supplied  angles. 

identities  the  type  of  aerosol  size  distribution  to  be  used, 
the  order  of  expansion  for  phase  functions  when  IANG  =  0, 
the  number  of  equally  spaced  angles  between  0°  and  180°  when 
IANG  =  1,  or  the  number  of  user  supplied  angles  when  IANG  =  2. 
used  for  optional  write/punch  on  NUNIT. 

the  back- scattering  coefficient  per  km  integrated  over  the 
size  distribution  for  each  WAVE. 

the  scattering  coefficient  per  km  integrated  over  the  size 
distribution,  for  each  WAVE. 
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SYMBOL 


Definition 


KTSUM 

NCRDS 


NINDX 

NUNIT 


OL(l) 

PFACT 

PSUM(J) 

SCAT(J) 

SFACT 

WAVE 


the  extinction  coefficient  per  km  integrated  over  the  size 
distribution,  for  each  WAVE. 

=  1  for  write/punch  only  Legendre  coefficients  on  NUNIT:  =  2 
write  only  phase  functions  and  scattering  fractions  on  NUNIT: 

=  3  write/punch  both  on  NUNIT. 

the  number  of  aerosol  components  which  will  have  different  optical 
constants,  mass  density  or  mass  concentration, 
defines  the  device  on  which  the  input  and/or  output  data  may  be 
stored  in  lieu  of  actual  card  punching;  may  be  used  to  place 
nominal  card  output  into  data  files  on  tape,  etc.  The  default 
value  (NUNIT  =0)  is  4. 

the  first  coefficient  in  the  Legendre  expansion  of  phase  functions. 
When  OL(l)  disagrees  with  ALBEDO  by  more  than  0.01  percent 
program  AGAUSX  should  be  rerun  using  a  larger  value  of  IT. 

=  WAVE*WAVE/(PI*CTSUM*EMM*EMM) .  When  PSUM(J)  is  multiplied  by 
PFACT  we  get  the  original  phase  functions. 

the  average  intensity  factors  integrated  over  the  size  distribution 
when  received  by  AGXPT3,  but  they  are  the  original  phase  functions 
for  each  WAVE  when  returned  to  MAIN. 

contains  'IT'  values  of  scattering  fractions  for  each  WAVE. 

=  WAVE*WAVE*DENS*10“6/4*PI*PI.  When  PSUM(J)  is  multiplied  by 
SFACT,  one  gets  the  scattering  fractions  of  ACT  code, 
the  wavelength  in  microns. 
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Subroutine  AGXPT3  -  Simplified  Flowchart 


RETURN 


1QQ 


SUBROUTINE  MIEGX* 


Subroutine  MIEGX  computes  various  efficiency  factors,  and  intensity  factors  i^ 
and  ±2  for  each  complex  refractive  index  m,  size  parameters  a  and  also  the 
first  three  Legendre  coefficients  (uQt  £5^  w2)  when  IAPX  is  greater  than 
zero.  The  Ricatti-Bessel  functions  and  their  derivatives  in  equations  (5)  and 
(6)  are  computed  by  the  forward  recursion  method.  The  initial  values  used  in 
forward  recursion  are: 


ii>0(z)  =  sin  z, 


^(z) 


sin  z 
z 


cos  z , 


XQ(z)  =  cos  z,  and 


X  _  cos_z  +  s^n  z 

-L  rr 


[Note:  £n(z)  =  ^(z)  +  i^Cz)] 


*A  Mie-type  routine  utilizing  a  continued  fractions  method  of  calculation  was 
obtained  from  W.  J.  Lentz  of  the  Atmospheric  Sciences  Laboratory  and  modified 
slightly  to  substitute  directly  for  MIEGX  in  AGAUSX,  thereby  eliminating  the 
need  for  the  forward  recursion  of  the  Ricatti-Bessel  functions.  The  new 
routine  is  more  accurate  for  larger  Mie  sizes  and  larger  imaginary  parts  of 
the  index  of  refraction.  This  new  routine  has  been  validated  by  New  Mexico 
State  University  under  contract  DAAD07-79-M-6275..  The  rotation  has  been  kept 
the  same  between  the  two  routines  as  much  as  possible. 


1Q1 


The  angular  functions  tt^  and  x^  are  also  computed  by  forward 
recursions  from  Equations  (24)  and  (25) .  The  initial  values  used  are 
7ro(0)  =  0,  7^(0)  =  1,  tq(0)  =  0,  and  x^E)  =  cos0. 

The  Mie  series  is  terminated  either  when  two  successive  terms  have 
(lRe(an)l  +  |lm(an)l  +  i Re (t*n)  |  +  |lm(bn)|)  <10  5,  or  when  the  number  of 
terms  exceeds  (8  +  Fa) .  F  is  1.2  for  a  £  51  and  is  1  +  2.26a-0'613  for 
a  >  51.  The  tolerance  value  of  10  3  could  be  decreased  for  higher  precision, 
if  the  user  so  desires. 

The  subroutine  computes  and  stores  arrays  of  a  and  b  until 

n  n 

convergence  is  reached  and  then  generates  necessary  and  x^  functions. 
Finally  it  computes  values  of  Q^,  Q^,  (absorption)  i^ 

(backscatter) ,  p(0),  and  radiation  pressure 


SYMBOL 

Explanation  or  Definition 

ALPHA 

Mie  size  parameter,  a  =  2irrA. 

C(K) 

the  array  of  cosine  of  the  scattering  angles.  There  are  'IT' 

elements  in  the  array. 

i 

CAY 

the  ratio  of  the  imaginary  part  to  the  real  part  of  adjusted 

refractive  index.  See  AG9PT2. 

EM 

the  real  part  of  adjusted  refractive  index. 

EN 

floating  point  representation  of  N. 

EYE1(K) 

ijL (Eq.  15)  at  angles  =  Arc  cos  [C (K) ] .  'IT'  elements  in  the 

array. 

EYE2 (K) 

i2(Eq.  16)  at  angles  =  arc  cos  [C (K) ] .  'IT'  elements  in  the 

array . 

FACT 

determines  the  cutoff  criterion  to  terminate  the  Mie  series. 

It 

is  equal  to  1.2  if  a  £  51  and  is  1  +  2.26a-0'613  for  a  >  51. 

Mie  theory  text. 

See 
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SYMBOL 


Definition 


FAN (N)  =  Im(a  );  Equation  (5). 

FBN(N)  =  Im(bn);  Equation  (6). 

GAMMA  the  true  imaginary  part  of  adjusted  refractive  index. 

IERROR  a  flag  which  is  set  to  unity  if  the  number  of  terms  (N)  reaches 
the  maximum  value  allowed  in  the  dimensions  of  a^  and  b^. 

ISW1  a  switch  parameter  used  in  applying  cutoff  criterion. 

01STAR  =  to^ ;  first  order  coefficient  for  Legendre  expansion  of  the 

average  intensity  P(K). 

02STAR  =  (^2')  second  order  coefficient  for  Legendre  expansion  of  the 
average  intensity  p(K). 

P(K)  =  (i  +  i  ) /2 ,  the  average  intensity  at  angles.  Arc  cos  (C(K)). 

1  2 

There  are  'IT'  elements  in  the  array. 

PIN  =  tt  (9);  Equation  (17). 

REAN  =  Re(an);  Equation  (5). 

REBN  =  RE(bn);  Equation  (6). 

RN  =  Re(^n(na));  Equation  (7). 

RNL1  =  ReC^fna));  Equation  (7). 

RPN  =  Re(^'(na)).  Prime  indicates  differentiation  with  respect  to 

the  argument. 

SGA  the  absorption  efficiency  factor. 

SGMAS  the  average  value  of  Cos (9),  where  9  is  the  scattering  angle. 
SGMP  the  radiation  pressure  factor  Qpr* 

SGR  the  back- scattering  efficiency  factor. 

SGS  the  scattering  efficiency  factor. 

SGT  the  extinction  efficiency  factor. 

SN  =  Im(^n(na));  Equation  (7). 
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SYMBOL 

SNL1 

SPN 

SUMI1I 

SUMI2I 

SUMI1R 

SUMI2R 

SUMRI 

SUMRR 

SUMS 

SUMS1 

SUMS2 

SUMT 


Explanation  or  Definition 


Im(^o(na));  Equation  (7). 


=  Im(Tjj'  (not) ) .  Prine  indicates  differentiation  with  respect  to 
the  argument. 

the  imaginary  part  of  the  scattering  amplitude  S^D); 

lmS  (0)  =  J  [Im(a  )7T  (0)  +  Im(b  )T  (0)]. 

1  L  n(n+l)  n  n  n'  n  ' J 

the  imaginary  part  of  the  scattering  amplitude  S2(0); 

ImSl(0)  ■  I  [I"(V,rn(6)  *  I"Can)T  (9)]. 


the  real  part  of  the  scattering  amplitude  S  (0)  ; 


Re  Sl(e)  ■  l  [R®tan)7rn(6)  *  Re(bn)Tn(e)]. 


the  real  part  of  the  scattering  amplitude  S2 (0) ; 

ReS2te)  ’  l  Sini  +  Re(an)Tn(0)]. 

N 

=  l  (-l)n(2n+l)Im(a  -b  ). 
n=l  n  n 

N 

=  l  01)n(an+l)Re(a  -b  ). 
n=l  n  n 

=  J (2n+l) ( | a  | 2  +  |b|2). 


-  I 

n=2 


(n-1) (n+1) 


[Re(an)Re(an+1)  +  Re(bn)Re(bn+1) 


+  Im(an)Im(an+1)  +  Im(bn)Im(bn+1)]. 


HTirfr  [Re(an)Re(bn)  +  Im(an)  ImCbJ  ] . 


=  l  (2n+l)Re(a  +  b  ) 
n  n 
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SYMBOL 


TAUN 

TERMN 

TN 

TNL1 

TPN 

UN 

UNL1 

UPN 


_ Explanation  or  Definition _ 

=  t  (0) ;  Equation  (17). 

=  |  Re  (an)  |  +  |lm(an)|  +  |Re(t>n)|  +  |  Im  (b^)  |  .  TERMN  determines 

the  cutoff  criterion  for  terminating  Mie  series. 

=  ^n(°0;  Equation  (7). 

=  (a) ;  Equation  (7). 

=  ^'(a);  prime  indicates  differentiation  with  respect  to  the 
argument . 

=  Xn(“);  Equation  (8). 

=  XQ (a) ;  Equation  (8). 

=  X^(°0;  prime  indicates  differentiation  with  respect  to  the 
argument . 
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RETURN 
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SUBROUTINE  ANGLE 


Subroutine  ANGLE  is  called  when  IANG  =  1  or  2.  It  is  usd  to  compute  the 
values  of  'IT'  equally  spaced  scattering  angles  between  0°  and  180°  or  a  read 
in  set  of  user-supplied  angles.  It  places  these  angle  values  in  the  array, 
H(I).  It  also  computes  the  cosines  of  those  angles  and  places  them  in  the 
array  C(I). 

SYMBOLS  _ Explanation  _ 


C(I)  the  array  containing  the  cosines  of  scattering  angles 

in  the  array  H(I). 

H(I )  the  array  containing  ! IT !  scattering  angles  (in  degrees). 

IT  is  the  number  of  scattering  angles  chosen  between  0°  and 

180° . 

RADS  =PI/180. 

SUBROUTINE  GUSET 


Subroutine  GUSET  uses  the  Davis  and  Rabinowitz  algorithm^  to  choose  n  values 
of  cos9kn  (k  =  l,2,...n)  between  the  interval  -1  <  cos9  <  1  and  the 

corresponding  values  of  quadrature  weights  akn.  The  abcissas,  cosg^  are  the 
n  zeros  of  the  Legendre  polynomials  Pn(cos9kn),  while  the  weights  are  given  by 


kn 


=  2(1  -  x 


kn 


(65) 


Up.  Davis  and  P.  Rabinowitz,  1956,  J  Res  NBS,  56:35 
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Initial  estimates  of  the  zeros  are  obtained  from  the  n  successive  zeros  of  the 
Bessel  function  (jQ(jk)  =  0)  via 


xkn  =  cos  +  j)2  +  (1  -  (~)  2)  /  4)  1/2] .  (66) 


final  values  of  the  xkn  are  found  by  Newton-Raphson  iteration.  The  tolernace 
of  Legendre  polynomial  zeros  is  set  at 


SYMBOL  Explanation 

AKN(k)  a^n  (Eq.  32).  'IT'  elements  in  the  array. 

n;  it  is  the  order  of  Legendre  expansion  for  phase  functions. 

P (N)  Legendre  polynomials,  Pn(x)« 

X  xkn  (Eq.  33). 

XKN(K)  cos6.  ;  p  (cos0.  )  =  0  within  10  ^ 
kn  n  kn 

TOL  the  tolerance  of  zeros  of  Legendre  polynomials  =  10 

Z(I)  jk  (Eq.  33).  I,  k  =  1,2,... IT. 

SUBROUTINE  GAUS 

Subroutine  GAUS  computes  the  Legendre  polynomials  as  PL(,)  and  computes 
Legendre  expansion  coefficients  given  by  equation  (31)  numerically.  To  do 
that  the  integral  in  equation  (31)  is  replaced  by  sunmation  as  follows: 


u>„  = 


_  (2fe+l) 


l  P(.e 


k- 1 


kn 


)p£(cosekn)akn. 


where  @kn  and  akn  have  the  meaning  as  given  in  subroutine  GUSET. 
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Using  the  values  of  coefficients  u>  and  equation  (30),  the  phase  functions  are 
reconstructed  and  called  Pc(9^n)*  GAUS  then  computes  the  root  mean  square 
deviation  between  the  original  phase  functions  p(0kn)  and  the  reconstructed 
phase  functions  Pc(9kn)  as  each  successive  term  is  added  to  the  series  in 
equation  (30).  Finally  it  prints  out  the  values  of  coefficients  and  rms 
deviations. 


SYMBOL 


F 

IT 

NCRDS 


NUN  IT 


OL(LL) 

PL(,) 

PC(I) 


RMS(J) 


Explanation  or  Definition 


is  the  array  containing  the  original  phase  function  P(0]tn). 

is  the  order  of  Legendre  expansion  for  phase  functions. 

=  1  for  write/punch  only  Legendre  coefficients  on  NUNIT:  =  2 
write  only  phase  functions  and  scattering  fractions  on  NUNIT: 
=  3  write/punch  both  on  NUNIT. 

defines  the  device  on  which  the  input  and/or  output  data  may 
be  stored  in  lieu  of  actual  card  punching;  may  be  used  to 
place  nominal  card  output  into  data  files  or  tape,  etc.  The 
default  value  (NUNIT=0)  is  4. 

the  array  containing  coefficients  wfc. 

Legendre  polynomials  (equivalenced  in  MIEGX). 

the  array  containing  the  reconstructed  phase  functions 
PC(0kn)* 

the  rms  deviation  between  the  original  phase  functions  and 
the  reconstructed  phase  functions. 


109 


SUBROUTINE  WATER 


The  purpose  of  subroutine  WATER  is  to  relieve  the  user  of  program  AGAUS  of  the 
task  of  looking  up  and  keypunching  data  on  the  index  of  refraction  and  mass 
density  of  liquid  water.  This  routine  receives  the  wavelength  (pm)  and 
temperature  (°K)  from  the  calling  program  as  variables  WVD  and  TEMPD.  It 
returns  the  mass  density  through  variable  DENSD,  and  the  real  and  imaginary 
parts  of  the  refractive  index  of  liquid  water  through  the  variables  EMD  and 
CAYD,  respectively. 

The  data  on  optical  constants  coded  into  routine  WATER  were  taken  from  the 
tabulation  by  Irvine  and  Pollack^  and  the  water  densities  were  taken  from  a 
copy  of  the  CRC  Handbook  of  Chemistry  and  Physics.  Tabulated  values  of  the 
real  (m)  and  imaginary  (k)  parts  of  the  refractive  index  are  available  for  the 
wavelength  range  0.20pm  to  200pm,  and  are  entered  at  uniform  wavelength 
increments.  Values  of  the  real  (m)  and  imaginary  (k)  parts  of  the  refractive 
index  at  wavelengths  other  than  found  in  the  table  are  estimated  through 
straight-line  interpolation.  Linear  interpolation  is  also  used  between 
tabulated  temperatures  in  calculating  the  mass  density  pw  (  =  DENSD). 

Methods  Used 


Subroutine  WATER  conducts  separate  binary  searches  of  the  wavelength 
table  LAMBDAQ  and  temperature  table  TEMPO  to  find  the  indices  L  and  L+l 
which  bracket  the  received  wavelength  (WVD)  and  temperature  (TEMPD).  It  then 
uses  linear  interpolation  to  get  estimated  values  of  EMT,  (M),  CAYT,  (k) ,  and 
RHODEN  (pw). 

The  interpolation  formula  can  be  written  as 


y(x)  =  y(x  )  +  [  At1  &](x-x  ) , 


with  y  =  m,  k  or  pw  and  x  =  wavelength  or  temperature. 


^W.  M.  Irvine  and  J.  B.  Pollack,  1968,  "Infrared  Optical  Properties  of  Water 
and  Ice  Spheres,"  Icarus ,  8:324 
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SYMBOL 

CAYD 

CAYT 

DENSD 

EMD 

EMT 

H 

L 

LAMBDA  ( 

NSUBI (  ) 

NSUBR 

P 

POINT 
TEMP (  ) 

TEMPD 

TMCHIJR 

RHODEN 

WAVE 


_ Explanation  or  Definition _ 

imaginary  part  of  refractive  index  k. 

the  interpolated  value  of  k  in  single  precision  form;  used 
intermediately  to  hold  summed  quantities. 

interpolated  result  for  the  mass  density  of  liquid  water  -  . 

real  part  of  index  of  refraction  -  m  . 

the  interpolated  result  for  the  real  part  of  the  refractive 
index  . 

an  indexing  (integer)  parameter, 
an  integer  indexing  parameter. 

array  of  wavelengths  at  which  data  are  entered  for  m  and  k; 
[typed  as  "real"]. 

array  of  values  for  k  (or,  nimaginary) 5  [typed  as  "real"], 
array  of  data  entries  for  m  (or,  nreai^  [typed  as  "real"], 
ah  (integer)  indexing  parameter, 
an  (integer)  indexing  parameter. 

array  of  temperature  values  (°K)  at  which  entries  for  pw 
exist  . 

temperature  at  which  pw  is  to  be  found  . 
temperature  at  which  value  of  pw  is  desired. 

interpolated  result  for  p^  , 

single  precision  version  of  wavelength  at  which  values 
of  n  and  k  are  desired. 
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SUBROUTINE  AGXPRT 


Subroutine  A3XPRT  prints  out  all  the  averaged  (sum/NWAVE)  results  essentially 
in  the  same  way  as  AGXPT3  does  for  each  wavelength  if  IAPX  =  0.  If  IANG  =  0, 
subroutine  GAUS  is  called  to  generate  Legendre  expansion  coefficients  which 
are  used  to  reconstruct  the  averaged  phase  functions.  GAUS  also  computes  the 
root  mean  square  deviation  between  the  original  phase  function  and  the 
reconstructed  phase  functions.  See  the  section  on  subroutine  GAUS  for  more 
details.  If  IANG  *  1  or  2,  the  above  calculation  is  skipped. 


SYMBOL  _ Explanation  or  Definition  _ 

C(I)  the  cosines  of  scattering  angles.  'IT'  elements  in  the  array, 

CATTN  the  average  (sum/NWAVE)  attenuation  coefficient  in  square  meters 
per  milligram  of  aerosol  material. 

H (I)  the  array  containing  'IT'  scattering  angles  (in  degrees). 

IANG  =  0,  for  the  computation  of  phase  functions  at  'IT'  Gauss- 

Legendre  quadrature  angles. 

=  1,  for  the  computation  of  phase  functions  and  scattering 
fractions  at  'IT'  equally  spaced  angles  between  0°  and  180°. 


IT 


KBAKT 


KEXTT 


KSCAT 


the  order  of  Legendre  expansion  for  phase  functions  when  IANG  = 

0,  or  the  number  of  equally  spaced  angles  between  0°  and  180° 
when  IANG  =  1,  or  the  number  of  user-supplied  angles  when  IANG  =  2. 
the  average  (sum/NWAVE)  back-scattering  coefficient  per  km, 
integrated  over  the  size  distribution, 

the  average  (sum/NWAVE)  extinction  coefficient  per  km,  integrated 
over  the  size  distribution. 

the  average  (sum/NWAVE)  scattering  coefficient  per  km, 
integrated  over  the  size  distribution. 
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SYMBOL 

NUNIT 

NWAVE 

PSUM(J) 

SCATT(J) 


_ Explanation  or  Definition _ 

defines  the  device  on  which  the  input  and/op  output  data  may  he 
stored  in  lieu  of  actual  card  punching;  may  be  used  to  place 
nominal  card  output  into  data  files  on  tape,  etc.  The  default 
value  (NUNIT=0)  is  4. 

the  number  of  wavelengths  or  relative  humidity  values  to  be 
treated  in  a  given  run. 

the  array  containing  the  values  of  average  (sum/ NWAVE) 
phase  function  integrated  over  the  size  distribution, 
the  array  containing  the  values  of  average  (sum/ NWAVE) 
scattering  fractions. 
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Subroutine  AGXPRT  -  Simplified  Flowchart 


RETURN 


114 


SUBROUTINE  GPHASX 


Subroutine  GPHASX  constructs  the  GHG  analytic  phase  function.  The  required 
input  consists  of  the  phase  function  at  zero  degrees,  oiQ,  u> 2>  and  the  order  of 
Legendre  expansion,  L.  The  parameter  g,  denoted  by  G1  in  the  subroutine,  is 
then  determined  by  an  iterative  process:  note  0  <_  g  <_  1.  We  then  determine  a 
and  subsequently  proceed  to  construct  the  GHG  phase  function,  checking  to 
insure  that  the  phase  function  does  not  go  negative;  if  it  does,  we  leave  the 
current  loop  and  starting  from  the  previous  value  of  p  (for  which  the  phase 
function  was  positive),  continue  to  construct  the  GHG  phase  function  from  the 
modified  HG  phase  function.  The  GHG  phase  function  is  then  renormalized  and 
the  appropriate  Legendre  coefficients  are  constructed. 


SYMBOL  _ Explanation  or  Definition _ 

ALPHA  .5[6J2/(&)oG2)-1]  . 

C  P(l)/w  . 

o 

COEF  Array  containing  first  three  Legendre  coefficients  (o)q,  o^,  w  )  . 

DELL  Iterative  increment  in  determination  of  G 

DEN  1  +  G2  -  2yG;  y  =  cosG. 

DENOM  (1  +  G2  -  2yG)1/2. 

ELHS  G [ (1  -  G2)P(1)/gjo  -  1  +  3G/2]  . 

EMO  Array  containing  Gauss-Legendre  angles. 

EMU  Array  containing  cosines  of  Gauss-Legendre  angles. 

ERR  Error  in  integrated  GHG  phase  function. 

G  Parameter  in  GHG  phase  function. 

GPFN  Array  containing  GHG  phase  function. 
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SYMBOL 


Explanation  or  Definition 


GSQ 

ITR 

NUN  IT 

OM 

OMO 

PA 

PHASEO 

PL 

PSUM 

RHS 

SUM 

TRUEG 

W 


G  *  G. 

number  of  iterations  in  the  determination  of  G;  100  iterations 
maximum. 

Auxiliary  unit  for  writing/punching  Legendre  coefficients  for 
GHG  phase  function. 

Array  containing  GHG  phase  function  Legendre  coefficients. 
Approximate  co  . 

Array  of  temporary  values  needed  in  subroutine  GUSSTX. 

GHG  phase  function  at  zero  degrees. 

Array  containing  Legendre  polynomials. 

Array  containing  true  phase  function  at  zero  degrees;  PSUM(l) . 

.5  w_/a)  . 

I  o 

integrated  GHG  phase  function. 

Asymmetry  factor  =  W-/(3a)  ). 

1  o 

Array  containing  weights  for  Gauss -Legendre  integration. 
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SUBROUTINE  GUSSTX 


This  subroutine  is  essentially  the  same  as  GUSET:  please  see  GUSET  for 

details. 

SUBROUTINE  PRINTX 


Subroutine  PRINTX  prints  out  various  quantities  for  the  analytic  GHG  phase 
function.  It  will  also  punch  out  the  Legendre  coefficients  for  the  GHG 
function. 


SYMBOL  _ Explanation  or  Definition  _ 

COEFS  Array  containing  GHG  Legendre  coefficients. 

EMU  Array  containing  cosine  of  Gauss -Legendre  angles. 

ERROR  Percent  error  in  the  GHG  integrated  phase  function. 

G  Asymmetry  factor;  0)^/(30)o). 

NUNIT  Auxiliary  unit  for  writing  (punching)  GHG  Legendre  coefficients. 
PHASE  GHG  phase  function  at  zero  degrees. 

PHFN  GHG  phase  function. 

PSUM  Array  containing  true  phase  function  at  zero  degrees. 

SUM  GHG  integrated  phase  function. 
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USERS  MANUAL  FOR  AGAUSX 


INPUT  DATA 

There  are  nominally  four  data  and  control  cards  required  for  running  program 

AGAUSX. 

CARD  _ General  Description _ 

I  A  set  of  integers  which  select  certain  options  available  within 
the  program. 

IA  A  user  set  of  supplied  angles. 

II  A  set  of  data  describing  the  parameters  of  the  particular  size 
distribution  function  to  be  used.  For  cases  other  than  the 
user-supplied  "arbitrary"  distribution  (IDSTP  =  0),  only  one 
card  of  this  type  is  needed. 

IIIA  A  set  of  data  describing  the  initial  wavelength  (ym)  to  be 

used,  its  increment,  aerosol  number  density  (particles  per  cc) , 
relative  humidity  (percent),  atmospheric  temperature  (degrees 
centigrade) ,  and  the  desired  convergence  testing  level. 

IIIB  If  the  looping  option  for  relative  humidity  is  invoked  then 
NWAVE-1  cards  containing  relative  humidity  and  atmospheric 
temperature  must  be  added. 

IV  A  set  of  data  describing  the  optical  and  physical  properties 

of  the  aerosol  material.  If  the  aerosol  is  a  mixture  of 
materials  of  unlike  properties  (NINDX  >  1),  more  than  one  card 
of  this  type  is  needed  for  each  card  of  type  III.  No  card  of 
this  type  is  used,  however,  for  runs  with  parameter  IDSTP 
equal  to  6. 
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Remarks : 


1.  The  simplest  type  of  run  (IDSTP=6,  water  cloud  or  fog  model) 
requires  one  card  each  of  types  I,  II,  and  III.  Forother  IDSTP  choices,  at 
least  one  card  of  type  IV  is  also  needed. 

2.  If  the  run  is  to  use  several  wavelengths  (NWAVE  >  1),  then  at 
least  one  card  of  type  IV  is  required  for  each  wavelength  (more  than  one  type 
IV  card  will  be  needed  at  each  wavelength  if  the  aerosol  is  a  multicomponent 
mixture  -  NINDX  >  1). 

3.  Two  special  "looping"  modes  of  AGAUSX  affect  the  number  of  cards 
of  type  III  which  are  needed : 

a.  For  runs  at  constant  relative  humidity  and  several 
wavelengths,  only  one  type  III  card  is  permitted.  To  activate  this  mode,  the 
wavelength  increment  DWAVE  on  card  type  III  must  be  larger  than  10“\ 

b.  Runs  at  constant  wavelength  and  a  set  of  differing  values  of 
relative  humidity  require  one  card  of  types  IIIA  and  IV  for  the  first  value  of 
relative  humidity;  subsequently  NWAVE-1  cards  of  types  IIIB  and  IV  are  needed 
for  each  additional  value  of  relative  humidity.  This  option  is  invoked  by 
setting  the  parameter  "DWAVE"  less  than  10-^  on  the  first  card  of  type  IIIA. 

Description  of  Types  of  Distributions 


IDSTP  DESCRIPTION 


0  This  is  an  arbitrary  user-supplied  distribution.  ' NRADI'  +  1 

cards  will  be  read;  the  first  card  contains  RLO  and  DELLR  (pm). 
The  rest  of  the  cards  carry  the  values  of  F(J)  and  must  be  in 
order  of  increasing  radius  value. 

1  This  is  the  zero-order  log-normal  distribution:  the 

distribution  function  is  given  by 
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F  (R)  = 


1 


/2t r  &oge(a)R 


iloge(R/R) 
&og e(a)  ^ 


1/2 


}. 


R  =  RBAR;  a  =  SIGMA,  is  the  standard  deviation.  This 
distribution  type  requires  one  input  data  card  to  read  in  the 
values  of  R,  a,  RLO  and  RHI. 

2  This  is  called  the  double  exponential  distribution  and  its 
distribution  function  is  given  by  F(R)  =  QA  exp(-AR)  + 

(l-Q)B  exp(-BR).  Q  =  CUE.  This  distribution  type  requires  one 
input  data  card  to  read  in  the  values  of  RLO,  RHI,  Q,  A,  and  B. 
Q  is  dimensionless  while  A  and  B  have  units  of  pm-1. 

3  This  model  (Deirmend jian' s  "Model  C")  does  not  require  any 
input  data  card.  It  carries  fixed  value  of  DENS,  RLO,  and 
DELRD.  RHI  is  determined  by  the  input  parameter  NRADI. 

F(R)  =  450.2  R  <  .08 

-  2.251*DELRD*R-4  R>_  .08 

4  The  distribution  function  of  this  model  (Junge  distribution) 
is  given  by  F(R)  =  QR-A,  Q  =  CUE.  This  distribution  type 
requires  one  input  data  card  to  read  in  RLO,  RHI,  A,  A. 

5  The  distribution  function  for  the  Modified  Gamma/ Generali zed 
Khirgian-Mazin  distribution  is 


F (R)  =  Ra  exp  [-(|-)Y  •  a  =  ALF,  y  =  GAM,  and 
c  ^ 
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Rc  =  RC.  One  input  card  is  needed  to  readin  RLO,  RHI,  RC, 

ALF,  GAM,  and  ELWC.  ELWC  is  not  needed  for  type  5  distribution 
and  therefore  can  be  left  blank. 

The  size  distribution  model  (NMSU  Fog  or  Cloud  Model)  is  very 
similar  to  type  5,  except  that  the  user  must  supply  one 
additional  input  parameter — namely,  the  liquid  water  content 
(ELWC)  in  gm/cc.  This  model  can  be  used  for  treating  situations 
involving  liquid  water  aerosols  like  clouds  or  fogs.  For  type 
6  runs  one  does  not  need  to  read  in  the  values  of  EMA,  CA¥A, 
RHOA,  CONC. 


This  distribution  is  essentially  same  as  Junge's  distribution 
(type  4)  except  that  it  has  fixed  parameters.  One  input  card 
is  needed  to  read  in  VIS  (visibility  in  km);  VIS  is  used  in 
calculating  DENS. 

This  is  a  fixed  parameter  Continental  Bi -modal  Model.  It  does 
not  require  an  input  data  card. 

This  is  a  fixed  parameter  Maritime  Bi -modal  Model.  It  does  not 
require  an  input  data  card. 

This  is  a  fixed  parameter  Urban  Bi -modal  Model.  It  does  not 
require  an  input  data  card. 


This  is  a  user-supplied  Bi -modal  Model.  This  requires  one 
input  card  to  read  in  FOA,  FOC,  SGA,  SGC,  RBARA,  RBARC.  Types 
8,  9,  10,  and  11  use  the  sum  of  two  log-normal  distributions: 


2  N. 

F(R)  =  l  - - -  exp 

i=l  /2tt  &oge(c^)R 


1 

2 


£oge(R/R.) 


] 


2 


=  FOA,  R^  =  RBARA,  =  SGA  with  similar  meaning  for  N^, 
R^,  and  Note  that  in  Type  8,  9,  10  the  values  of  SGA 
and  SGC  are  £oge(a). 
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12 


This  model  (Marshall-Palmer  Bain  Model)  is  a  simple  exponential 
model  which  assumes  an  empirical  relationship  between  rain  rate 
and  droplet  size  distribution  parameters  : 


F(D)  =  NQexp(-AD). 


Nq  =  0.08  cm-^,  and  A  =  41(RN)-0*21  in  which  RN  =  RAIN  is  the 
rain  rate  in  mm/hour.  Diameter  D  is  in  cm.  The  corresponding 
size  distribution  function  of  radius  R  is  given  by 


F(R)  =  2NQexp(-2AR) . 


This  distribution  requires  one  input  data  card  to  read  in  RAIN. 
The  values  of  RLO  and  RHI  are  fixed  at  0  and  0.5  cm 
respectively.  Due  to  the  limitations  on  the  range  of  Mie- 
sizes  (subroutine  MIEGX)  type  12  usage  is  limited  to  wave¬ 
lengths  of  the  order  of  1  mm  or  larger.  Since  subroutine 
WATER  does  not  contain  optical  data  for  wavelengths  longer 
than  0.2  mm,  type  12  runs  require  the  user  to  supply  the 
values  of  EMA,  CA¥A,  and  RHOA  as  if  rain  were  a  non-aqueous 
aerosol. 
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SUMMARY  OF  DATA  CARD  REQUIREMENTS 


CARD 

TYPE 

I 


_ INPUT  SYMBOLS _ 

Integer  control  parameters:  NWAVE,  NINDX,  IW,  IDSTP, 

NRADI,  IT,  IANG,  NCRDS,  ITOT,  MQRTE ,  IAPX,  FORMAT  (1215). 

NWAVE:  is  the  number  of  wavelengths,  or  relative  humidity 
values  to  be  treated  in  this  run.  See  comments  circa  read 
of  WAVE,  DWAVE,  etc. 

NINDX:  is  the  number  of  aerosol  components  which  will  have 

different  optical  constants,  mass  densities  or  mass  concentrations. 

IW:  =0  will  set  the  refractive  index  of  the  dry  aerosol  equal 
to  that  of  water  at  the  input  wavelength  and  temperature  -  see 
Card  4,  EMA,  CAYA. 

IDSTP:  Identifies  type  of  aerosol  size  distribution  to  be  used. 

NRADI:  Number  of  particle  radii  to  be  expected  for  IDSTP  =  0  or 
3:  the  input  value  of  NRADI  is  ignored  for  IDSTP  not  zero  or  3. 
NRADI  must  be  less  than  513. 

IT:  is  the  number  of  Gauss -Legendre  angles  (order  of  expansion) 
if  IANG  =  0  or  number  of  angles  between  0  degrees  and  180  degrees 
if  IANG  =1.  If  only  extinction  coefficients,  etc.  are  desired, 
i.e.,  not  phase  functions,  then  set  -IT-  equal  to  one. 

IANG:  =  0  for  computations  of  phase  function  at  -IT-  Gauss- 
Legendre  quadrature  angles:  no  scattering  fractions  will  be 
printed.  IANG  =  1  for  computations  of  scattering  fractions  and 
phase  functions  at  -IT-  equally  spaced  angles  between  0  and 
180  degrees.  IANG  =  2  will  allow  -IT-  user  supplied  angles  to 
be  read  -  FORMAT  (16F5.1).  This  requires  at  least  one  card  of 
type  IA. 
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CARD 

TYPE 


INPUT  SYMBOLS 


NCRDS :  =  1  write/punch  Legendre  coefficients  and  order  of 

expansion  on  NUNIT  -  FORMAT  (D25. 14, IX, 15) :  =  2  write/punch 

scattering  fractions,  phase  functions,  cosines  and  counter 
on  NUNIT  -  FORMAT  (3(D12. 61X) ,15) :  =  3  write/punch  both 

scattering  fractions,  phase  functions,  cosines,  counter,  and 
Legendre  coefficients  and  order  of  expansion  on  NUNIT  - 
FORMATS  (3 (D12 . 6 . , IX) , 15)  and  (D25 . 14 , IX, 15) :  scattering 
fractions,  phase  functions,  cosines  and  counter  (respectively) 
are  always  written/punched  before  the  Legendre  coefficients 
and  order  of  expansion  (respectively) . 

ITOT:  Works  in  conjunction  with  NCRDS  and  NWAVE  to  write/punch 
either  individual  wavelength  values  or  averaged  wavelength 
values  or  both  on  NUNIT:  =  1  for  individual  wavelengths; 

=  2  for  averaged  wavelengths;  =  3  for  both. 

NUNIT:  defines  device  number.  See  previous  comment  card. 

MQRTE :  =  12345  will  cause  prints  of  Mie  efficiency  factors  at 

every  value  of  particle  radius  used  in  the  Mie  calculations; 
set  MQRTE  =  0  if  such  prints  are  not  desired. 

IAPX:  if  GT  zero  the  GHG  (analytic)  phase  function  will  be 

constructed  at  IAPX  Gauss-Legendre  angles  (order  of  expansion) . 
User  supplied  set  of  -IT-  angles  FORMAT  (16F5.1)  16  values  per 
card,  more  than  1  card  may  be  needed.  This  card  is  only  needed 
when  IANG  =  2. 
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CARD 

TYPE 

II 


_ INPUT  SYMBOL _ 

Distribution  parameters:  only  one  type  per  run.  FORMAT  (6E12.6). 
For  a  more  detailed  description  of  these  parameters  see  the  previous 
discussions  in  AGXPT1.  All  units  are  in  micrometers. 

TYPE  0.  USER  SUPPLIED  -  NRADI  +  1  CARDS 
RLO,  DELLR 
FF (I) ,  I  =  1,  NRADI 
TYPE  1.  LOG-NORMAL 

RBAR,  SIGMA,  RLO,  RHI 
TYPE  2.  DOUBLE  EXPONENTIAL 
RLO,  RHI,  CUE,  A,  B 
TYPE  3.  DEIRMENDJIAN  MODEL  C 

-  NO  INPUT  - 

TYPE  4.  POWER  LAW  (JUNGE) 

RLO,  RHI,  CUE,  A 
TYPE  5.  MODIFIED  GAMMA 

RLO,  RHI,  RC,  ALF,  GAM,  ELWC 
TYPE  6.  MODIFIED  GAMMA  FOG  MODEL 

RLO,  RHI,  RC,  ALG,  GAM,  ELWC 
TYPE  7.  POWER  LAW 
VIS 

TYPE  8.  CONTINENTAL  BIMODAL 

-  NO  INPUT  - 

TYPE  9.  MARITIME  BIMODAL 

-  NO  INPUT  - 
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CARD 

TYPE  _ INPUT  SYMBOLS _ _ 

II  TYPE  10.  URBAN  BIMODAL 

-  NO  INPUT  - 

TYPE  11.  USER  SUPPLIED  BIMODAL 

FOA,  RBARA,  GSA,  FOC,  RBARC ,  SGC 
TYPE  12.  MARSHALL- PALMER  RAIN  MODEL 
RAIN 

III  Control  Parameters:  FORMAT  (6E12.6)  WAVE,  DWAVE,  RELHUM, 

DENSH,  TEMP,  DELTA.  For  looping  over  different  values  of 
relative  humidity  set  DWAVE  =  0,  and  add  NWAVE  -  1  cards 
containing  RELHUM  and  TEMP.  FORMAT  (2E12.6). 

WAVE:  is  wavelength  in  micrometers 

DWAVE:  is  the  wavelength  increment  in  micrometers.  If 
DWAVE  is  less  than  l.E-4,  a  special  case  applies  used  for 
looping  over  NWAVE  values  of  RELHUM.  For  example  in  that 
case  a  card  which  carries  the  value  of  RELHUM  and  TEMP  must 
be  repeated  NWAVE  -  1  times;  the  first  card  of  this  type 
must  contain  WAVE,  DWAVE,  RELHUM,  DENSH,  TEMP  and  DELTA. 
RELHUM:  is  relative  humidity  in  percent. 

DENSH:  is  particle  number  per  cubic  centimeter.  User-supplied 
value  of  DENSH  will  be  ignored  for  IDSTP  =  3  or  greater  than 
6  because  those  distributions  carry  pre-determined  density 
values.  Also,  if  DENSH  is  less  than  IE-4,  the  particle  number 
density  will  be  calculated  from  mass  densities  and  mass 
concentrations . 
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CARD 

TYPE 


INPUT  SYMBOLS 


III  TEMP:  is  the  temperature  of  the  atmosphere  in  degree^  C. 

DELTA:  is  the  convergence  criterion.  Within  a  particular 

size  range  interval,  halving  is  terminated  when  the  quantity 
DEL  is  less  than  DELTA. 

IV  Optical  and  physical  data  (Read  in  AGXPT2)  FORMAT  (4F10.6E15.7) 
EMA,  CAY A,  EMUA,  RHOA,  CONC. 

EMA:  is  the  real  part  of  the  index  of  refraction  of  dry  aerosol. 
CAYA:  is  the  imaginary  part  of  refractive  index  for  dry  aerosol. 
CAYA  is  assumed  to  be  negative:  do  not  enter  CAYA  with  a 
negative  sign. 

EMUA:  is  Hanel's  growth  factor  (ME-BAR) /accretion  coefficient. 

RHOA:  is  the  mass  density  (sp.  GRAV)  of  dry  aerosol. 

CONC:  is  the  mass  concentration  (GM/CC)  of  dry  aerosol. 

For  looping  over  wavelength  or  relative  humidity  repeat  this 

card  NWAVE  times  and  interleave  with  card  3:  For  looping  over 

aerosol  components  also  repeat  nindx  times.  This  card  is  not 

needed  when  IDSTP  =  6.  Also  if  IW  =  0  EMA  and  CAYA  will 

internally  be  set  equal  to  the  refractive  index  of  water  and 

therefore  may  be  left  blank :  values  for  EMUA,  RHOA,  and  CONC 
Remarks’  maY  necessarY  to  suPPly>  depending  upon  the  users  needs. 

1.  Card  Type  IV  is  not  required  and  is  never  read-in  if  IDSTP  =  6 
(water  cloud/for  model)  since  the  relevant  data  are  obtained  from  subroutine 
WATER. 

2.  In  the  case  IDSTP  =  12,  the  rain  model,  card  type  IV  must  carry 
the  optical  data  for  liquid  water  as  EMA  and  CAYA.  The  reason  for  that 
inconsistency  is  that  the  rain  model  will  most  likely  be  useable  at 
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wavelengths  for  which  no  data  appear  in  subroutine  WATER  (restricted  to 
X  <  0.2mm).  [The  large  droplet  sizes  to  be  found  in  the  rain  model  will 
usually  cause  premature  failure  of  the  Mie  routine  at  small  wavelengths.] 

3.  Card  type  IV,  with  appropriately  adjusted  data  values,  must  be 
repeated  (NWAVE  -  1)  times. 

4.  EMA  and  CAYA  may  be  left  blank  on  Card  type  4  when  IW  =  0  as  they 
will  internally  be  set  equal  to  the  value  for  water. 

Incidental  Remarks  Regarding  AGAUSX 

(1)  In  AGAUSX,  the  average  "dry  volume  per  particle"  is  found  using 
all  available  values  for  particle  radii.  It  may  differ  from  that  printed 
in  routine  AGXPT2  since  that  routine  might  not  proceed  to  the  use  of  all 
available  values. 

(2)  The  "volume"  convergence  tests  in  AGXPT2  use  the  volume  inferred 
after  any  growth  arising  from  non-zero  saturation  ratios  has  been 
included. 

(3)  It  should  be  noted  that  the  convergence  tests  used  in  AGAUSX 
do  not  really  provide  tests  of  the  absolute  accuracies  achieved  for  the 
tested  quantity.  It  is  possible  for  "exit"  to  occur  0$<A)  even  though  the 
use  of  "another"  halving  might  lead  to  $>A  again.  If  runs  using  some 
choice  of  A  might  lead  to  what  appear  to  be  "unusual"  or  unexpected 
results,  it  is  advised  that,  as  a  test,  the  case  be  re-run  using  a  smaller 
value  of  A. 

(4)  Some  users  may  wish  to  explore  possible  increases  in  computation 
efficiency  which  might  result  from  changes  in  the  number  (NI)  of  size- 
intervals,  or,  the  ways  in  which  they  have  been  chosen  (in  AGXPT1) . 
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Comments  on  Usage  of  AGAUSX  in  the  IANG  =  0  Mode 


AGAUSX  has  been  coded  to  preserve  the  option  of  creating  the  expansion 
coefficients,  &)^,  which  are  needed  by  the  multiple  scattering  codes  STAR04 , 

AG SCAT ,  and  the  thermal  emission  code  CLEM70(8).  Usage  of  the  on  coefficients 

X/ 

for  the  above  purposes,  however,  requires  some  caution  on  the  part  of  the 
user.  In  particular,  users  must  not  assume  that  the  absence  of  abnormal 
termination  of  runs  of  AGAUSX  guarantees  that  "all  is  well".  Past  experience 
has  shown  that  "warnings"  of  possible  problems  printed  by  those  programs 
may  not  be  noticed  or  taken  very  seriously.  In  an  attempt  to  overcome 
those  oversights,  additional  tests  and  warnings  have  been  incorporated  in 
AGAUSX.  If  those  warnings  are  to  be  meaningful,  users  should  examine 
printed  outputs  carefully  for  such  warnings  before  using  the  results  of 
runs  as  inputs  to  subsequent  codes. 

One  particular  point  which  should  be  checked  is  to  see  that  the 
quantity  printed  under  the  heading  ALBDO  agrees  reasonably  well  with  the 
zeroth  coefficient  (o)q)  of  the  expansion  coefficients.  Substantial 
disagreement  between  those  two  quantities  usually  means  that  the  parameter 
"IT"  used  in  a  run  was  too  small  to  achieve  a  really  accurate  reconstruction 
of  the  phase  function  from  only  "IT"  Legendre  expansion  terms. 

Should  possible  inaccuracies  be  found,  please  contact  R.  C.  Shirkey, 
Atmospheric  Sciences  Laboratory,  White  Sands  Missile  Range,  New  Mexico 
88002  . 
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SAMPLES  OF  INPUT  AND  OUTPUT  FOR  AGAUSX 

INPUT 

The  subsequent  output  section  is  based  upon  the  following  input. 
Example  1 :  IDSTP  =  0  (User  supplied  distribution  data) 

Column  5 

2190  11  5000  012345  0 

0. 10000E-00  0.0200E+00 
0. 00000E-00 
0. 20000E-00 
0. 40000E-00 
0. 60000E-00 
0. 80000E-00 
1 . 00000E-00 
0. 80000E-00 
0. 60000E-00 
0. 40000E-00 
0. 200Q0E-00 
0. 00000E-00 

10. 6000E+00  0.10000E-05  00.0000E+00  0.00000E-05  25.0000E+00  1.00000E-02 

1.9530  .46800  0.00000  1.87000  1.000000E-09 

75. 0000E+00  25. 0000E+00 

1.9530  .468000  0.15900  1.87000  1.000000E-09 


13Q 


This  example  is  set  up  to  read  11  distribution  data  cards,  and  to  run  for  two 
values  of  relative  humidity  at  X  =  10.6pm.  Since  MQRTE  =  12345,  detailed  Mie 
results  for  each  of  the  11  radii  will  be  printed.  Parameter  'IT'  is  5, 

IANG  =  1.  Note  that  the  second  data  card  contains  RLO  (=  0.1pm)  and  DELLR 
(  =  0.02).  The  next  eleven  cards  contain  F(RLO) ,  F(RL0  +  DELLR,  etc.).  The 
remaining  cards  consist  of  two  pairs  of  cards,  two  for  each  of  the  two  values 
of  NWAVE.  In  all  cases,  DENSH  is  zero  which  means  that  particle  number 
densities  will  be  calculated  from  RHOA  (1.87,  above),  CONC  (10-9  gm/cc,  above) 
and  the  average  particle  volume  computed  within  the  program. 

Example  2:  IDSTP  =  1  (Log  normal  distribution) 


Column  5 

i 

139545000000 

3. 70000E-01  1 . 54000E+00  0.00500E-00  1.00000E-00 

1 0.  600QE+00  0. 10000E-05  00.000GE+00  O.OOOOOE-05  25.0000E+00  1.00000E-02 

1.9530  .468000  0.00000  1.87000  1.000000E-09 

0.15900  1.87000  1.000000E-09 

0.15600  1.87000  1.000000E-09 

Remarks  :  The  second  card  carries  RBAR,  SIGMA,  RLO  and  RHI.  If  RLO  and  RHI 
are  set  to  zero,  the  program  itself  will  choose  RLO  and  RHI. 
Particle  number  density  will  also  be  calculated  from  mass  density 
and  concentration. 
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Example  3  :  IDSTP  =  3  (Deirmend jian' s  Model  C) 


Column  5 

2103  100  5000000 

10. 6000E+00  0. 10000E-05  00.0000E+00  0.00000E-05  25.0000E+00  1.00000E-02 

0.0000  .000000  0.00000  1.87000  1.0000000-09 

75.00000+00  25.00000+00 

0.0000  .000000  0.15900  1.87000  1.0000000-09 

Remarks  :  100  values  of  radii  will  be  used  for  2  values  of  relative  humidity 

with  different  growth  factors.  Particle  density  will  be  calculated 
from  mass  density  and  concentration.  Since  IW  =  0,  the  aerosol  com¬ 
ponent  is  considered  to  be  pure  water;  i.e.,  EMA  and  CAYA  will  be 
set  equal  to  the  values  for  water. 

Example  4:  IDSTP  =  5  (Modified  gamma  distribution) 

Column  5 

~r 

1  2*95  40  1000000 

.0278520+00  1.225490+01  4.000000+00  .6000000+01  .1000000+01 

. 700000CH-00  .0000000-04  .0000000+00  1.000000+02  5.000000+00  1.000000-02 

1.3300  .000000  0.00000  1.00000  1.0000000+00 

1.4300  .000000  0.00000  1.00000  1.0000000-01 

Remarks :  This  example  sets  NINDX  =  Z,  and  treats  a  mixture  of  two  aerosol 
components  having  different  refractive  indices  and  mass  con¬ 
centrations. 
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Example  5:  IDSTP  =  7  (Power  Law) 

Column  5 

110701131  23  08 

+6.000000+00 

10.60000+00  .0000000+00  .0000000+00  .1000000+03  25.00000+00  1,000000-02 

0.0000  .000000  0.00000  0.00000  0.0000000+00 

Remarks :  The  GHG  analytic  phase  function  will  be  constructed  at  eight 

Gauss-Legendre  angles  (IAPX  =  8)  for  pure  water  (IW  =  0).  Scattering 
fractions  (not  computed),  phase  functions,  cosines,  a  counter, 
and  Legendre  coefficients  will  be  written  on  unit  23  for 
individual  wavelengths  only. 
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RELATIVE  HUHIDITT  FOR  THIS  RUN  ■  7S.00  PERCENT.  WAVELENGTH  ■  10.600  NJCRONS 

INDEX  OF  REFRACTION  FOR  PURE  WATER  IS?" I. 178200  .  .07S620I 
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integer  control  parameters:  nnave  ninox  iw  idstp  nradi  it  iang  ncros  itot  nunit  mqrte  iapx 
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type  (MICRONS)  INDEX  ( S«  MICRONS)  (SQ  MICRONS) 
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PROGRAM  AGAUSX 
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ro»  component  no.  2  interval  no.  i  33  radii  were  used,  contribution  to  ctsum  »  2.668959+001 

for  Component  no.  2  interval  no.  2  »7  radii  were  used,  contribution  to  ctsOm  -  1 .928298+002 

For  component  no.  z  :  vpf  -  6.17099-012  mass  concentration  ■  6.17099-012  sm/cc*  kext  ■  1.69517-003  per  km 


WARNING  OPTIMAL  PF  EXPANSION  ORDER  OF  M2  EXCEEDS  INPUT  IT  *  1.  PF  VALUES  SHOULD  BE  USED  CAUTIOUSLY 
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integer  control  Parameters:  N*aVE  nindx  iw  idstp  nradi  it  iang  ncrds  ITOT  NUNIT  MflRTE  I APX 
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CALL  GAUS ( PSUM , NUN ! T , I  TOT )  PT3X0660 
ROUTINE  GAUS  GENERATES  and  P R I N T S / P UNCHE S  THE  LEGENDRE  PT3X0670 
expansion  coefs  iomegasi  ior  the  phase  function.  '  PT3X0680 
check  TO  SEE  if  SNf-,.  SCAT.  ALBEDO  (  ALBDO  I  COMPUTED  DIRECTLY  PT3X0690 
FROM  CROSS-SECTIONS  AGREES  WITH  THAT  FOUND  FROM  THE  LEGENDRE  PT3X0700 
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